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ABSTRACT 

A  prerequisite  for  the  design  of  low  frequency  radio  receivers  is 
a  model  for  low  frequency  atmospheric  radio  noise  that  encompasse  the 
non- Gaussian  nature  of  the  actual  noise  process  and  is  sufficiently  tract¬ 
able  to  enable  performance  analysis  and  optimization  of  receiver  designs. 
This  work  describes  a  new  model  for  atmospheric  noise  waveforms 
observed  at  the  output  of  the  antenna  bandlimiting  filter.  This  model, 
which  is  based  on  statistical  analysis  of  sample  records  of  these  wave¬ 
forms,  is  used  to  analyze  the  performance  of  typical  radio  navigation 
receivers  and  to  determine  near  optimum  receiver  performance.  The 
analysis  is  verified  by  simulating  the  receiver  structure  and  testing 
the  receiver  with  the  actual  noise  sample  records. 

The  proposed  model  for  bandlimited  atmospheric  noise  utilizes  a 
background  Gaussian  process,  of  constant  power  level,  to  which  is  added 
discrete  bursts  of  Gaussian  noise  whose  power  level,  for  each  burst,  is 
a  random  variable.  The  burst  occurrences  and  the  power  level  are  com¬ 
mon  to  waveforms  observed  in  disjoint  frequency  channels  and  provides 
a  model  for  the  statistical  dependence  between  such  waveforms.  The 
model  is  shown  to  match  first  order  characteristics  of  a  bandlimited 
noise  waveform  and  a  noise  envelope,  recorded  in  disjo  at  channels, 
over  a  range  of  noise  conditions  from  nearly  Gaussian  to  he  severest 
noise  caused  by  frontal  thunderstorm  activity.  The  length  of  each  noise 
burst,  and  occurrence  rate  intensity,  are  both  stochastic  in  the  model. 

A  non-homogeneous  Markov  process  provides  a  description  of  these  burst 
characteristics  and  is  shown  to  result  in  an  envelope  autocorrelation  that 
matches  observed  correlations  over  periods  of  1  millisecond  to  1  second. 

THESIS  SUPERVISOR:  Amar  G.  Bose 
TITLE:  Professor  of  Electrical  Engineering 
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The  atmospheric  noise  model  is  used  to  analyze  the  performance  of 
three  basic  types  of  phase -locked -loops  used  in  navigation  receivers. 

The  first  of  these  is  a  linear  time -invariant  loop  optimized  for  the  noise 
power  level.  The  second  adds  a  zero  memory  nonlinearity  in  the  path  of 
the  loop  error  signal  and  is  optimized  to  the  first-order  noise  probability 
density.  T.ie  third  is  a  near  optimum  structure  that  uses  the  joint  chan¬ 
nel  noise  model  and  the  second  channel  noise  envelope  to  convert  the 
atmospheric  noise  process  to  an  approximate  time-varying  gaussian  pro¬ 
cess.  This  can  be  optimally  filtered  by  a  time -varying  Kalman-Bucy 
estimator.  It  is  shown,  using  the  recorded  noise  sequences,  that  the 
simple  nonlinear  loop  provides  performance  improvements  of  up  to  16  dB 
compared  to  the  time -invariant  linear  case,  and  that  the  optimum  time- 
varying  filter  provides  additional  improvements  of  up  to  3-5  dB  when 
receiver  processing  intervals  are  less  than  several  seconds. 

The  work  concludes  with  suggestions  for  further  research,  including 
an  application  of  the  joint  channel  noise  model  to  the  design  of  an  optimum 
digital  communication  receiver  that  appears  to  be  both  practical  and  offer 
significant  performance  advantages  over  present  designs. 
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The  Very-Low  Frequency  (VLF,  3-30  kHz)  and  Low-Frequency  (LF, 
30-300  kHz)  electromagnetic  spectrum  is  extensively  used  for  digital 
radio  communications  and  radio  navigation.  When  the  performance  of 
these  radio  systems  is  limited  by  noise,  such  noise  is  generally  atmo¬ 
spheric  radio  noise,  additively  combined  with  the  signal  at  the  receiving 
antenna.  Thermal  noise  in  the  antenna-receiver  system  is  generally  of 
a  lower  level  than  this  radio  noise.  The  principal  source  of  the  domi¬ 
nant  atmospheric  noise  is  radiation  from  lightning  discharges,  which 
propagates  considerable  distances  at  these  frequencies.  When  atmo¬ 
spheric  noise  waveforms  are  observed  at  the  output  of  the  receiver  band- 
limiting  filters  used  at  these  frequencies,  it  is  described  with  words 
such  as  "impulsive,"  "peaky,"  "crashy,"  etc.,  recognizing  the  fact 
that  individual  lightning  discharge  events  are  discernible  and  tend  to 
dominate  the  waveform  at  any  given  time  instant.  A  result  of  this  indi¬ 
vidual  dominance  is  that  statistical  descrip  tions  of  the  noise  waveform 
are,  in  general,  distinctly  non-gaussian  since  the  individual  dominance 
condition  clearly  violates  the  basic  assumption  of  the  central  limit  theo¬ 
rem.  The  design  of  receiver  structures  for  processing  signals,  cor¬ 
rupted  by  non-gaussian  atmospheric  noise,  is  a  difficult  procedure  and 
the  results  depend  critically  upon  the  statistical  characteristics  ascribed 
to  the  noise.  Improved  characterization/design  procedures  are  useful 
goals  for  VLF/LF  systems  because  of  the  large  costs  associated  with 
obtaining  equivalent  receiver  performance  improvements  by  increasing 
radiated  signal  power.  For  example,  the  costs  of  erecting  antenna 
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structures  vary  as  the  fourth  power  of  the  height  above  800-1000'. 

The  research  reported  here  has  addressed  the  problem  of  improving 
the  statistical  characterization  of  low-frequency  atmospheric  noise  for 
the  purpose  of  improving  the  design  of  VLF/LF  radio  systems.  The 
succeeding  work  is  divided  into  four  parts.  Chapter  1,  Introduction, 
provides  a  qualitative  description  of  the  physical  noise  process  to  moti¬ 
vate  our  conceptual  approach  to  the  problem.  This  approach  is  placed 
in  perspective  to  the  results  reported  by  other  workers  during  the  past 
twenty  years.  Chapter  2,  Experimental  Data,  is  a  tabulation  of  the  vari¬ 
ous  measurements  we  made  of  atmospheric  noise  waveforms  using  a 
digital  data  collection  system.  Most  of  this  data  represents  new  con¬ 
tributions  to  the  literature  of  this  field.  Chapter  3,  An  Atmospheric 
Noise  Model,  presents  a  mathematical  model  for  bandlimited  atmospheric 
noise  waveforms  which  describes  first-order  statistics,  time  structure, 
and  the  statistical  dependence  of  waveforms  observed  in  different  fre¬ 
quency  channels.  Chapter  4,  An  Optixnum  Design  for  Navigation  Re¬ 
ceivers,  applies  classical  servo-mechanism  theory  to  the  design  of  a 
nonlinear  Phase-Locked  Loop  (PLL),  and  exploits  the  co-frequency- 
channel  dependence  of  the  noise  to  design  a  near-optimum  Kalman-Bucy, 
time-varying  PLL.  Performance  predictions  of  the  designs,  based  on 
the  noise  model,  are  verified  with  the  use  of  atmospheric  noise  sample 
records.  Chapter  5  concludes  with  suggestions  for  further  work,  in 
particular,  with  other  applications  of  the  noise  model  to  digital  com¬ 
munication  receiver  design. 


1.  1  Atmospheric  Noise  Mechanism 
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Electromagnetic  radiation  from  lightning  discharges,  the  principal 
source  of  atmospheric  noise,  have  been  extensively  studied.  While 
there  is  a  wide  diversity  of  specific  details  in  these  measurements, 
there  is  agreement  as  to  the  general  structure  of  these  discharges.  A 
basic  cloud -to-ground  discharge  consists  of  a  series  of  short  current 
pulses  associated  with  the  advance  of  the  ionization  of  the  air  dielectric 
along  the  incipient  discharge  path,  followed  by  cloud -to-ground  move¬ 
ment  of  the  electric  charge  which  created  the  initial  potential  difference. 
Figure  1-1  (A)  shows  the  rauiation  from  such  a  basic  discharge  with  the 
ionization  noise  and  large  cloud -to -ground  pulse  clearly  observable. 

This  oscilloscope  photograph  was  made  at  the  amplified  output  of  a 

3 

broadband  loop  antenna  with  -3  dB  bandwidth  of  15  to  85  kHz.  Maxwell 
reports  that  the  spectral  energy  of  the  initial  ionization  burst,  termed 
the  leader  stroke,  "S  centered  at  30  kHz,  whereas  the  main  stroke, 
termed  the  return  stroke,  has  a  peak  of  spectral  energy  in  the  region 
of  3-8  kHz.  The  leader  stroke  has  been  observed  to  have  an  average 
length  of  1  msec,  whereas  the  return  stroke  falls  in  a  100-200  psec 
category. 

A  large  number  of  discharge  events,  perhaps  a  majority  under 
conditions  of  large,  high  energy  storm  cells,  are  composed  of  compli¬ 
cated  repetitions  of  the  basic  discharge  just  described.  Repeated  return 
strokes  are  termed  multiple  discharges  and  have  been  observed  to  con- 

i  4 

sist  of  20-30  such  strokes  lasting  for  periods  of  200-500  msec.  ' 

Figure  1-1  (B)  shows  several  return  strokes  occurr  ing  m  a  short  time 
frame,  while  1-1  (C)  shows  a  longer  segment  of  a  multiple  discharge. 
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Of  particular  interest  to  our  work  is  the  fact  that  these  more  complicated 
discharge  patterns  have  a  "leader"  type  pattern  that  becomes  a  continu¬ 
um  noise  burst  rather  than  a  series  of  identifiable  strokes  or  events. 
Since  this  leader  structure  has  a  principal  spectral  energy  above  30  kHz, 
whereas  the  return  strokes  are  centered  below  this  frequency, we  would 
expect  significant  differences  in  the  time  structure  of  noise  waveforms 
observed  at  different  frequencies  with  respect  to  30  kHz.  This  is  indeed 
often  the  case,  as  is  seen  in  Figure  l-l(D).  This  tendency  of  complex 

discharge  events  to  appear  as  continuous  noise  bursts  has  been  noted  by  * 

4  5  4 

many  experimental  workers.  *  Gupta  provides  comparative  tracings 
of  single  and  multiple  discharge  noise  patterns  observed  at  the  output  of 
a  narrow-band  filter  followed  by  an  envelope  detector.  These  tracings, 
made  at  frequencies  in  the  MF  and  HF  region,  clearly  show  the  compli¬ 
cated  burst  structure  of  both  the  single  and  multiple  discharge  phenom¬ 
ena. 

We  have  been,  up  to  this  point,  implicitly  considering  the  structure 
of  lightning  discharges  observed  near  the  source  where  we  may  assume 
single  mode  (ground-wave)  propagation  from  source  to  observer.  In 
determining  the  total  observed  atmospheric  noise  spectrum  one  must  also 
consider  propagation  effects  which  are  very  significant  and  change  radi¬ 
cally  from  VLF  to  LF,  and  above.  The  sum  total  of  these  effects  is 

3 

given  by  Maxwell  in  determining  the  broad  spectral  characteristics  of 
atmospheric  noise.  All  of  the  work  in  this  area  has  been  summarized 
by  Oh  who  provides  a  plot  from  kHz  to  GHz  of  relative  noise  intensity 
and  the  references  for  measurements  at  various  frequencies.  We  shall 
not  be  concerned  with  these  relative  levels  of  atmospheric  noise;  rather. 
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we  have  sought  to  provide  a  qualitative  background  on  the  nature  of  the 
discharge  process  to  motivate  our  conceptual  approach  to  the  problem 
and  to  point  out  the  fundamental  changes  in  structure  that  occur  with 
frequency.  These  changes  and  the  effect  of  propagation  will  be  seen  in 
many  of  the  measurements  to  be  reported  in  Chapter  2. 

1.2  Experimental  Noise  Data 

The  statistical  measurement  of  atmospheric  radio  noise  that  has 
been  most  often  reported  in  the  literature  is  the  exceedance  probability 
(one  minus  the  cumulative  probability)  of  the  noise  envelope  at  the  output 
of  a  bandlimiting  filter.  A  representative  but  not  exhaustive  list  of  such 
work  is  given  by  references  1,2,5,7-13,30.  The  major  reason  for 
the  use  of  this  characterization  is  its  ease  of  measurement  with  rela- 

O 

tively  simple  analog  circuits,  threshold  detectors,  and  counters.  Evans 
and  Coon  recorded  noise  waveforms  on  analog  tape  recorders  and  later 
derived  these  measurements  from  digital  computer  processing  of  the 
sample  records.  The  exceedance  probability  is  usually  plotted  on  Log- 
Log  versus  Log  paper  where  a  Rayleigh  distribution  appears  as  a 
straight  line  of  slope  2.  From  these  measurements  workers  have  con¬ 
cluded  that  the  atmospheric  noise  waveform  consists  of  a  low  amplitude 
Rayleigh  component  (implying  a  low -amplitude,  gaussian-distributed 
component  in  the  narrow -band  waveform)  with  a  large  amplitude  com¬ 
ponent  due  to  distinct  or  dominating  lightning  bursts.  Most  workers  have 
used  a  log-normal  or  l/y°  probability  density  function  to  fit  the  observed 
large  amplitude  excursions. 

The  only  statistical  characterization  of  the  time  structure  of 
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atmospheric  noise  waveforms  that  has  been  reported  has  been  some 
type  of  level  crossing  measurement.  Generally  this  takes  the  form  of 
a  probability  density  for  the  time  until  a  given  level  is  crossed  in  the 
upward  direction  by  the  noise  envelope,  given  that  an  upward  crossing 
occurred  at  time  zero.  The  principal  conclusion  from  these  measure¬ 
ments  has  been  that  the  discrete  discharge  events,  assumed  to  be  de¬ 
tected  by  such  a  measurement,  are  not  independently  distributed  in  time 
(a  Poisson  process)  but  rather  tend  to  cluster,  i.  e. ,  there  is  a  higher 
than  expected  probability  that  another  event  will  occur  after  an  initial 
event  has  been  detected.  This  effect  is  a  manifestation  of  multiple  dis¬ 
charge  phenomena. 

In  considering  what  type  of  statistical  parameters  to  measure  to  char 
acterize  the  time  structure  of  atmospheric  noise,  one  must  bear  in  mind 
that  a  theoretically  infinite  set  of  probability  distributions  is  required  to 
completely  describe  this  structure.  The  usefulness  of  any  particular 
parameter  selected  can  only  be  judged  in  relation  to  its  end  use.  While 
level  crossing  measurements  are  quite  illuminating  of  the  temporal  be¬ 
havior  of  discharge  events,  it  is  not  clear  that  they  are  equally  useful 
to  the  receiver  engineer,  especially  at  higher  band-center  frequencies 
where  the  complex  leader  structure  becomes  dominant.  As  will  be  seen 
in  Chapter  2,  we  have  selected  a  much  different  characterization  of 
atmospheric  noise  time  structure. 

1 . 3  Proposed  Noise  Models 
I.  3.  1  Impulse  Model 

A  number  of  different  models  for  atmospheric  noise  waveforms  have 
appeared  in  the  literature  over  the  past  15  years.  The  majority  of  these 
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might  be  termed  "impulse  models"  for  they  derive  their  motivation 
from  an  approximation  of  the  actual  generating  mechanism  by  an  impulse 
train  exciting  the  narrow-band  observation  filter.  One  might  also  term 
them  physical  models  because  the  model  parameters  can  usually  be 
closely  linked  to  physical  quantities  or  filter  parameters.  The  classic 
model  of  this  type  is  that  given  by  Furutsu  and  Ishida.***  Their  model 
includes  not  only  a  Poisson  distributed  wave  of  impulse  functions  of  ran¬ 
dom  energy,  but  further  Poisson  distributed  pulse  bursts  following  a 
basic  impulse.  They  were  concerned  with  MF  and  HF  atmospheric  noise 

and  this  type  of  model  accounts  for  the  complex  leader  or  continuum 

10  14  15 

structure  dominant  at  these  frequencies.  The  general  approach  '  ’ 

to  the  determination  of  noise  statistics  at  the  output  of  the  observation 
filter,  for  this  type  of  model,  has  been  to  use  Price's  technique  with 
characteristic  functions  to  determine  the  first  order  and  higher  proba¬ 
bility  densities  in  the  characteristic  function  domain.  The  resulting 
expression  cannot  in  general  be  transformed  to  yield  a  probability  den¬ 
sity  function  (pdf)  except  at  the  large  and  small  amplitude  limits.  The 

2 

noise  model  given  by  Beach  is  of  this  impulse  type  and  is  the  most  com¬ 
plex  model  that  we  are  aware  of,  including  a  number  of  different  rate 
parameters  for  weather  patterns,  time-of-year,  etc.  Beckman** 
utilizes  a  different  approach  with  much  the  same  result  in  terms  of  the 
envelope  pdf.  These  models  are  well  suited  for  the  inference  of  physical 
parameters  (such  as  storm  cell  distribution  or  propagation  parameters) 
from  observed  noise  waveforms  but  are  difficult  to  use  for  receiver 
engineering,  nor  do  they  provide  insight  into  the  actual  noise  structure 
on  which  the  receiver  processor  must  operate.  An  exception  to  this 


seems  to  occur  for  ELF  noise  waveforms  as  considered  by  Evans  and 

Modestino  where,  due  to  the  frequency  (3-300  Hz)  and  bandwidth,  the 

noise  is  completely  dominated  by  the  return  stroke  phenomena  which  at 

1 6 

this  bandwidth  is  well  modeled  by  an  impulse  function.  Snyder  also 
uses  an  impulse  model  as  a  vehicle  for  the  derivation  of  the  optimum 
estimator  equations  for  a  binary  detection  receiver.  Snyder's  approach 

r 

is  only  dependent  upon  structure,  avoiding  any  explicit  reference  to 

probability  distributions.  However,  the  resulting  equations  form  an 

infinite  set  of  coupled  partial -differential-integral  equations  which  are 

1 5 

as  difficult  to  exploit  as  the  determination  of  the  complete  model 
statistics  at  the  filter  output. 

1.3.2  Communication  Models 

A  number  of  atmospheric  noise  models  have  appeared  in  the  com- 
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munications  literature  and  have  been  directed  towards  receiver 

design.  A  principal  feature  of  these  models  is  that  they  attempt  to 

describe  the  noise  waveform  at  the  output  of  the  bandlimiting  filter  with 

only  peripheral  concern  for  the  underlying  generating  mechanism. 

1 7 

Kapp  utilized  a  "switched  process"  model,  originally  due  to  Kurz, 

in  which  the  observed  waveform  was  given  as  either  a  bar.dlimited 

gaussian  process  or  an  impulsive  process  with  first  order  Cauchy  pdf. 

The  probability  of  "finding"  one  or  the  other  of  these  distributions  is 

then  a  measure  of  the  impulsiveness  of  the  noise  waveform. 

1 7 

A  much  different  model  was  given  by  Shaft  in  which  the  envelope 
of  the  narrow-band  process  is  given  as  EXP(n(t)),  where  n(t)  is  a 
gaussian  process.  This  expression  only  provides  correct  large  amplitude 


-25- 

behavior  of  the  model;  however.  Shaft  reports  satisfactory  application 
of  the  model  to  the  prediction  of  VLF  modem  performance. 

1.3.3  Multiplicative  Noise  Model 
18 

Hall  proposed  a  model  in  which  the  narrow-band  noise  waveform 
at  the  filter  output  is  given  as  y(t)  =  A(t)  •  n(t),  where  n(t)  is  a  conven¬ 
tional  gaussian  process  and  A(t)  is  a  lowpass  random  process.  The  time 
structure  of  n(t)  is  that  obtained  by  illuminating  the  bandlimiting  filter 
with  white  noise.  The  modulating  process,  A(t),  is  used  to  increase  the 
dynamic  range  of  the  basic  gaussian  process  to  match  that  exhibited 
by  atmospheric  noise  waveforms.  At  the  same  time,  the  time  structure 
of  A(t)  is  available  to  imitate  the  complex  time  behavior  of  the  dis¬ 
charge  process.  Hall  demonstrated  that  with  an  assumed  generating 
mechansim  (which  will  be  covered  in  detail  in  Chapter  3)  the  first  order 
statistics  of  either  the  envelope  or  the  narrow-band  rf  process  could 
be  expressed  in  closed  form  for  his  model  and  that  this  form  matched 
quite  well  the  various  reported  exceedance  probability  distributions  of 
the  noise  envelope. 

Recalling  the  discussion  of  the  complex  nature  of  the  discharge 
event,  in  particular,  the  tendency  at  LF  and  above  for  the  noise  mech¬ 
anism  to  be  dominated  by  the  leader  or  fine  structure,  we  see  that 
there  is  considerable  motivation  for  the  multiplicative  noise  model. 

The  modulating  process,  A(t),  we  can  associate  physically  with  the 
fluctuating  power  level  of  the  atmospheric  noise  at  a  receiving  site, 
caused  both  by  different  sources  and  the  long  noise  burst  characteris¬ 
tics  of  multiple  discharges  at  individual  sources.  The  use  of  A(t)  to 
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to  model  these  effects  at  the  filter  output  offers  the  possibility  of  approx¬ 
imating  the  very  complex  behavior  of  the  lightning  discharge  process  by  a 
simpler  description  of  a  random  generating  mechanism  for  A(t). 

Another  concept  that  we  have  evolved  m  considering  this  type  of 
model  is  that  A{t),  modeling  the  envelope  of  the  broadband  discharge 
process,  should  be  common  to  noise  waveforms  observed  in  nearby  but 
disjoint  frequency  channels.  Noting  that  our  final  goal  is  the  use  of  the 
model  in  receiver  design,  this  concept  has  considerable  importance  for 
if  one  knew  A(t)  exactly,  then  the  corrupting  noise  waveform  would  be¬ 
come  a  time-varying  or  non-stationary  gaussian  process.  Existing  sig¬ 
nal  processing  designs  generally  include  time-varying  gaussian  noise 
as  the  most  general  case,  and  hence  the  difficult,  non-gaussian,  atmo¬ 
spheric  noise  problem  would  have  the  potential  for  being  converted  to  a 

m"'Li  more  tractable  time -varying  gaussian  problem.  While  there  have 
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been  a  few  references  in  the  literature  *  to  using  information  from 
adjacent  frequency  channels,  these  are  based  on  an  assumed  determinis¬ 
tic  relationship  between  noise  waveforms,  that  both  are  generated  by 
an  ideal  impulse  exciting  the  respective  filters.  Our  postulated  model 
appears  to  be  the  first  attempt  at  modeling  a  statistical  relation  between 
such  waveforms. 

We  began  the  experimental  phase  of  this  research  by  approaching 
the  design  of  a  data  collection  system  to  answer  certain  questions  about 
the  suitability  of  such  a  conceptual  multiplicative  model.  The  specific 
measurements,  the  system  designed  to  obtain  these,  and  the  analyzed 
results  are  given  in  the  next  chapter. 


Chapter  2 

EXPERIMENTAL  DATA 
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In  Chapter  1  we  discussed  the  nature  of  electromagnetic  radiation 
from  lightning  discharges  as  it  affects  low-frequency  radio  noise.  From 
this  discussion  and  a  review  of  past  efforts  at  modeling  the  noise  pro- 

1 8 

cess,  we  postulated  a  multiplicative  noise  model  as  suggested  by  Hall. 
This  model  is  given  as  y{t)  =  A(t)  •  n(t),  where  y(t)  is  the  output  of  a 
bandlimiting  filter  excited  by  received  atmospheric  noise  fields  and 
represents  the  waveforms  that  must  be  processed  by  the  receiver.  The 
term  n(t)  is  a  conventional  narrow-band  gaussian  process  and  A(t)  is  a 
lowpass  random  process.  This  conceptual  form  of  model  had  two 
appealing  features  for  our  requirement  of  a  noise  model  that  was  par¬ 
ticularly  useful  for  receiver  design:  1)  The  modulating  process  A(t) 
could  be  used  to  absorb  the  very  complex,  long-time  nature  of  the  actual 
discharge  process  which  is  broadband  and  excites  all  nearby  frequency 
channels;  2)  a  knowledge  of  or  good  approximation  of  A(t)  would  convert 
an  untractable  non-gaussian  noise  problem  into  a  much  simpler  time- 
varying  gaussian  problem  for  which  optimum  solutions  are  known.  In 
this  chapter  we  describe  an  experimental  program  formulated  to  answer 
the  two  basic  questions  raised  by  point  1)  above:  Are  the  actual  noise 
waveforms  amenable  to  these  interpretations?  Our  purpose  in  this 
chapter  is  to  provide  a  broad  overview  of  the  measured  noise  charac¬ 
teristics  prior  to  detailed  development  of  a  mathematical  model  and 
comparison  of  model  characteristics  with  the  observations. 


2.  1  Noise  Characteristics  Measured 
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The  instrumentation  system  that  was  constructed  for  this  work  is 
described  in  section  2.  2  below.  This  system  was  designed  to  measure 
the  following  statistical  characteristics  of  bandlimited  atmospheric  noise 
waveforms  at  frequencies  of  14,  65,  and  83  kHz  and  bandwidths  of  1,  10, 
and  20  kHz. 

1)  Probability  Densicy.Noise  waveforms  were  sampled  at  periods  of 

5 

1  to  20  milliseconds  with  sample  records  of  10  samples.  These  were 
analyzed  to  estimate  the  noise  probability  density  (pd)  using  a  rectangular 
approximation  over  small  intervals. 

2)  Joint  Probability  Density.Simultaneous  samples  of  a  noise  wave¬ 
form  in  frequency  channel  one  and  the  noise  envelope  in  a  non -overlapping 
frequency  channel  two  were  recorded.  These  were  analyzed  to  estimate 
the  joint  pd  surface  and  various  conditional  statistics  of  these  two,  depen¬ 
dent,  variables. 

3)  Conditional  Sampling.Samples  of  noise  waveforms  were  taken, 
conditioned  by  hardware  upon  a  given  event  in  another  channel,  to 
reinforce  the  measurements  of  (2). 

4)  Autocorrelation.  The  autocorrelation  of  the  noise  waveform  was 
measured  by  sampling  the  waveform  at  an  integral  multiple  of  the  band- 
center  frequency  period  and  autocorrelating  the  resulting  record. 

5)  A(t)  Estimates.  To  generate  an  estimate  of  the  short-term  power 
level  of  the  atmospheric  noise,  represented  by  A(t)  in  the  model,  we 
integrated  the  absolute  value  of  the  noise  waveform 
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for  fixed  time  periods  of  1  and  10  msec.  It  is  well  known  that  a 
statistically  optimum  estimate  of  an  unknown,  gaussian  noise,  power 
level  is  given  by  an  average  of  the  observations  squared.  However,  for 
atmospheric  noise  waveforms  with  a  dynamic  range  in  excess  of  80  dB 
this  requires  squaring  circuitry  with  a  dynamic  range  of  160  dB,  which 
is  presently  unattainable.  The  use  of  the  magnitude  was  selected  as  a 
realizable  alternative.  These  estimates  of  A(t^)  were  made  in  two  fre¬ 
quency  channels  simultaneously  and  used  to  determine  the  auto  and 
cross-channel  correlation  of  A(t). 

The  exact  computational  algorithms  used  to  realize  these  estimates 
and  approximate  error  analyses  are  given  in  Appendix  B. 

2.  2  Instrumentation  System 

A  functional  description  and  performance  measurement  of  the  in¬ 
strumentation  system  is  given  in  Appendix  A.  The  basic  functions  of 
the  system  are  def,  mbed  below. 

1)  Analog  Processing.  Atmospheric  noise  fields  were  received  on 
a  one-meter  loop  antenna,  amplified,  and  filtered  with  a  bank  of  LC 
notch  filters  to  remove  possibly  saturating  communication  signals.  The 
notch  filtered  signal  was  bandlimited  to  four  channels  by  fixed  tuned  LC 
filters.  The  dynamic  range  of  the  system  to  this  point  was  greater  than 
100  dB,  with  system  noise  better  than  6  dB  below  the  quiescent  or  back¬ 
ground  atmospheric  noise  level.  These  four  channels  were  processed 
through  two  wide -band  amplifier  samplers  or  two  reclify-integrate-hold 
circuits.  The  sampler  dynamic  range  was  greater  than  80  dB,  with  a 
600 -nsec  gate, and  included  appropriate  offset  adjustments.  The  active 
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rectifier  circuits  had  a  linear  range  of  76  dB.  The  rectifier  outputs 
could  also  be  fed  to  a  20-kHz  lowpass  filter  and  thence  to  the  wide-band 
samplers  for  envelope  data  collection. 

2)  Digital  Processing.  Digital  processing  began  with  a  1  MHz  ultra¬ 
stable  crystal  oscillator,  which  was  used  as  a  system  clock  to  control 
all  sampling  and  data  handling  functions.  A  wide  variety  of  sampling 
modes  could  be  generated,  from  rates  of  100  kHz  to  0. 1  Hz  and  in  con¬ 
tinuous  or  burst  patterns.  Special  circuits  allowed  coarse  synchroniza¬ 
tion  of  the  sampling  pattern  with  received  pulse  signals  from  the  Loran-C 
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navigation  system  so  that  noise  data  could  be  taken  at  this  operating 
frequency  (100  kHz)  without  signal  contamination.  The  four  sampled 
analog  channels  could  be  connected  to  the  Analog-to-Digital  converter, 
in  various  periodic  sequences,  to  provide  15  bit,  ±10  volt  digitizing 
capability  to  all  channels.  Digitized  data  was  stored  in  an  8K  by  16  bit 
core  memory  prior  to  asynchronous  recording  on  paper  tape. 

3)  Data  Recording.  Digital  data  was  recorded  on  a  IK  bit-per-second 
paper  tape  punch  for  convenient  input  to  the  small  computer  system  used 
throughout  this  research  program.  Data  flow  was  asynchronous  from 

the  buffer  memory  via  an  optical  isolator  to  prevent  transient  contami¬ 
nation  of  the  analog  system. 

A  complete  description  of  parts  (2)  and  (3)  of  the  data  system  is 
24 

given  by  Lee. 

2.  3  Classification  of  Observations  by  Weather  Conditions 

The  instrumentation  system  was  extensively  tested  on  known  deter¬ 
ministic  and  random  signal  sources  as  described  in  Appendix  C.  Prelim- 
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inary  recording  and  analysis  of  atmospheric  noise  was  conducted  at 
Cambridge,  Massachusetts  during  July  1971.  The  complete  data  record 
on  which  the  remainder  of  this  work  is  based  was  made  during  the  period 
10  August  to  22  August  1971  at  the  United  States  Coast  Guard  Electronics 
Engineering  Center,  Wildwood,  New  Jersey,  which  is  located  90  miles 
from  the  nearest  urban  area.  The  loop  antenna  was  mounted  on  top  of 
a  single  story  concrete  block  building  approximately  l/4  mile  from  the 
ocean.  The  antenna  reception  null  zones  were  oriented  to  minimize 
received  signals  from  the  Loran-C  transmitting  stations  at  Nantucket, 
Massachusetts,  and  Cape  Fear,  North  Carolina. 

During  the  data  collection  period  the  Eastern  Seaboard  experienced 
three  distinct  types  of  weather  patterns  that  were  highly  correlated  with 
the  observed  atmospheric  noise  characteristics.  For  convenience  in 
referring  to  the  different  noise  conditions,  we  have  termed  the  three 
basic  weather/noise  conditions  "quiet,"  "tropical,"  and  "frontal."  The 
frontal  condition  occurred  on  11  August  1971  and  was  caused  by  an  ad¬ 
vancing  cold  front  with  extensive  thunderstorms  along  the  squall  line. 
Heavy  rain  and  visual  lightning  passed  over  the  observation  site  3  hours 
after  the  frontal  data  was  recorded.  The  weather  during  the  next  week 
was  dominated  by  a  high  pressure  air  mass  that  had  caused  the  frontal 
conditions.  This  air  was  cool  and  stable  with  very  little  thunderstorm 
activity  and  produced  the  quiet  conditions.  As  this  weather  pattern 
moved  to  the  east,  a  llow  of  warm,  moist,  unstable  air  from  the  Gulf 
of  Mexico  overspread  the  East  Coast.  This  tropical  weather  period  had 
a  high  geographic  density  of  thunderstorm  activity,  especially  in  the 
afternoon  at  the  time  of  maximum  accumulated  energy  in  the  vertical 
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air  thermals.  As  will  be  seen,  the  quiet  noise  conditions  were  quite 
close  to  being  gaussian  and  hence  these,  or  gaussian  conditions,  repre¬ 
sent  probably  a  majority  of  year  'round  noise  conditions  temperate 
latitudes.  The  tropical  conditions  we  would  describe  as  typical  "dog 
days  of  August"  weather  that  one  experiences  during  the  East  Coast 
summer.  The  frontal  conditions,  producing  severest  noise  conditions, 
were  only  observed  at  five  distinct  times  during  the  July-August  period 
that  we  made  observations  and  were  all  associated  with  frontal  weather 
patterns.  Frontal  noise,  then,  represents  atypical  temperate  latitude 

behavior,  but  is  possibly  representative  of  more  common  conditions  in 

1'! 

equatorial  regions.  ’  In  addition  to  these  three  basic  noise  conditions, 
a  transitional  condition  was  ob,-  »,veo  luring  the  night  when  the  tropical 
air  flow  was  moving  over  the  n-aldie  Atlantic  region.  This  is  labeled 
"quiet-night"  in  what  follows. 

Our  sample  of  atmospheric  noise  conditions,  while  taken  in  a  rela¬ 
tively  small  time  frame,  is  probably  representative  of  a  wide  variation 
of  noise  conditions  that  one  would  find  throughout  the  world.  This  state¬ 
ment  will  be  given  more  meaning  in  Chapter  3  where  we  shall  observe 
that  the  various  noise  conditions  can  be  related  to  extreme  ranges  in 
important  model  parameters. 

2.  4  First  Order  Probability  Density  of  Atmospheric  Noise 

2.4.  1  Single  Frequency  Channel  Probability  Density 

The  range  of  the  first  order  pd  of  a  bandlimited  atmospheric  noise 
waveform,  at  VLF  and  LF,  is  shown  in  Figure  2-1.  These  observations 
were  made  at  each  frequency  with  different  gains  which  attempted  to 
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optimize  use  of  the  80  dB  dynamic  range  of  the  A-D  converter.  The 
combined  plots  were  made  by  removing  the  differential  scale  factor, 
resulting  in  a  pd  that  would  have  been  measured  by  a  fixed  gain  system  of 
greater  dynamic  range.  The  scale  factor  between  frequency  channels 
is  arbitrary,  no  attempt  was  made  in  this  work  to  determine  relative 
noise  levels  between  frequency  bands  (see  Maxwell,0  Oh  ). 

The  most  striking  factor  in  these  plots  is  the  wider  variation  in 

pd's  at  LF  (65  kHz)  than  at  VLF  (14  kHz).  This  is  a  consequence  of  two 
3 

physical  factors,  propagation  and  the  frequency  dependence  of  the  dis¬ 
charge  characteristics.  Propagation  attenuation  is  a  factor  of  10  greater 
(@)  8000  KM)  at  65  kHz  compared  to  14  kHz,  increasing  to  a  factor  of  100 
at  100  kHz.  Thus  the  received  atmospheric  noise  field  at  VLF  is  sen¬ 
sitive  to  a  much  wider  geographic  area  than  at  LF.  At  the  New  Jersey 
observation  site  this  meant  that  the  14  kHz  observations,  during  quiet 
conditions,  included  discharge  effects  from  much  farther  south,  regions 
of  greater  thunderstorm  activity  than  affected  the  65  and  100  kHz  obser¬ 
vations.  During  the  frontal  observations,  the  received  noise  field  was 
dominated  by  the  intense  lightning  activity  along  the  cold  front  squall 
line.  At  this  time,  the  long  multiple  discharge  phenomena  produced 
a  near  continuum  of  noise  at  LF  and  above,  while  the  VLF  noise  con¬ 
sisted  principally  of  distinct  pulses  from  the  return  stroke  (see  Fig¬ 
ure  1  -  1(D)) .  These  observations  support  our  qualitative  analysis  of 

expected  differences  in  VLF  and  LF  noise  waveforms. 

6  18  25 

Several  workers  have  suggested  ’  *  that  a  bandwidth  of  10  kHz 

plays  an  important  roll  in  the  behavior  of  atmospheric  waveforms.  The 
arguments  given  involve  the  dependency  of  peak  value  of  a  noise  bu?st 
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on  bandwidth,  with  10  kHz  being  the  point  at  which  this  dependence 
changes  from  a  linear  relation  l.o  a  square  root  relationship.  This 
prediction  is  based  on  a  comparison  of  the  average  length  of  a  return 
stroke  (100  jisec)  to  the  effective  filter  impulse  response  time,  and  a 
consideration  as  to  whether  the  discharge  appears  as  a  noise  burst  or 
an  impulse.  It  is  not  clear  what  effect  this  consideration  has  on  the 
noise  pdf,  however,  our  observations  at  LF  show  no  important  differ¬ 
ence  in  noise  pd's  measured  at  1,  10,  and  20  kHz  bandwidths.  Fig¬ 
ure  2-2  shows  the  noise  pd  at  65  kHz,  in  1  and  10  kHz  bandwidths, 
plotted  on  a  normalized  amplitude  scale  illustrating  this  fact. 


Figure  2-2.  LF  Noise  pd's  Observed  in  1  and  10  kHz  Bandwidths 
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We  believe  this  behavior  is  attributable  to  the  importance  of  the 
leader  stroke  component  of  the  noise  burst,  which  for  a  single  discharge 
has  an  average  length  of  1  msec.  One  would  suspect,  if  this  is  the  case, 
that  observations  made  in  bandwidths  much  less  than  1  kHz  would  show 
significant  differences,  however,  we  have  not  made  any  observations 
at  these  bandwidths.  We  shall  discuss  the  role  of  filter  bandwidths  less 
than  1  kHz  in  Chapter  3  as  they  significantly  affect  signal  processing 
structures.  However,  in  the  remainder  of  this  chapter,  we  shall  treat 
all  of  the  LF  data  as  equivalent  at  the  three  observations  bandwidths  of 
1,  10,  and  20  kHz  with  only  notations  on  the  figures  to  indicate  the  band¬ 
width  used.  All  of  the  data  recorded  at  14  kHz  was  taken  in  a  1  kHz 
bandwidth. 

2, 4.  2  Joint  Frequency  Channel  Probability  Density 

As  we  have  noted,  a  principal  feature  of  the  multiplicative  noise 
model  that  we  have  postulated  is  that  the  modulating  random  process, 
A(t),  is  common  to  nearby  but  disjoint  frequency  channels.  To  explore 
this  concept  we  recorded  samples  of  the  bandpass  noise  in  one  frequency 
channel,  termed  the  rf  channel,  and  the  bandpass  envelope  in  an  adjacent 
frequency  channel,  termed  the  envelope  or  pilot  channel.  The  decision 
to  use  the  envelope  as  a  short  term  measure  of  the  fluctuating  power 
level  rather  than  the  A(t^),  used  to  estimate  the  dynamic  properties  of 
A(t),  was  based  on  tests  to  determine  which  of  these  parameters  had  the 
highest  correlation  coefficient  with  the  magnitude  of  the  bandpass  noise. 
The  results  of  these  tests  are  shown  in  Figure  2-3  where  the  use  of  the 
envelope  is  seen  to  result  in  a  slightly  larger  cross -correlation  coeffi¬ 


cient. 
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Figure  2-3.  Correlation  of  RF  Noise  Magnitude  with  Various  A(t)  Esti¬ 
mates 
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Figure  2-4  shows  two  oscilloscope  photographs  of  the  pilot  channel 
envelope  and  rf  channel  noise  for  a  multiple  discharge  and  scattered 
low  level  single  discharges.  Inspection  of  the  relationship  of  the  pilot 
channel  envelope  (dots  were  added  to  distinguish  this  trace)  shows  that 
this  envelope  does  not  "follow"  the  rf  noise  level  in  a  simple  deter¬ 
ministic  manner.  For  example,  the  four  large  rf  bursts  labeled  A, 

B,  C,  and  D  do  not  have  the  same  response  in  the  pilot  channel.  In 
principle,  if  we  knew  the  state  of  the  respective  rf,  and  the  pilot  chan¬ 
nel  bandlimiting  filters  at  some  time  and  the  exact  form  of  the  broadband 
atmospheric  noise  excitation  for  future  time,  we  could  calculate  these 
responses  deterministically.  The  implicit  assumption  of  our  multipli¬ 
cative  noise  model  is  that  we  do  not  have  this  information;  we  seek 
rather  a  noise  model  that  describes  these  relationships  in  an  average 
or  statistical  sense.  It  would  appear  from  these  photographs  that  our 
concept  of  a  multiplicative  noise  model,  with  A(t)  independent  of  fre¬ 
quency,  would  provide  a  plausible  basis  for  modeling  these  relations. 

RF  channel  and  pilot  channel  waveforms,  such  as  those  seen  in 
Figure  2-4,  were  sampled  at  1  kHz  rates  in  bursts  of  approximately  10 
seconds.  The  burst  length  was  determined  by  input/output  rate  and  mem¬ 
ory  capacity  of  the  instrumentation  system.  The  pilot  channel  sample 
was  taken  0.7  msec  after  the  rf  sample  to  compensate  for  the  differing 
group  delays  in  the  respective  bandlimiting  filters.  This  delay  differ¬ 
ence  produced  the  lag  in  the  cross-correlation  peak  seen  in  Figure  2-3 
and  the  obvious  delay  seen  in  the  photographs  of  Figure  2-4.  The  pilot 
channel  center  frequency  was  83  kHz  with  a  0.9  kHz  bandwidth  for  all 
recordings,  whereas  rf  channels  were  located  at  65  kHz  with  1  and 
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and  10  kHz  bandwidths  and  100  kHz  with  a  20  kHz  bandwidth.  Appen¬ 
dix  A  provides  a  description  of  these  filters.  For  the  100  kHz  record, 

sampling  bursts  were  restricted  to  30  msec  in  length  and  were  time- 
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interlaced  with  the  on-air,  pulsed,  Loran-C  navigation  signals.  Joint 
channel  sample  records  were  not  made  at  VLF  frequencies. 

These  joint  channel  sample  records  were  used  to  estimate  the  joint 
probability  density  surface  for  rf  samples  from  channel  1  and  envelope 
samples  from  channel  2.  The  analysis  algorithm  is  given  in  Appendix  B. 
The  probability  density  surface  is  shown  in  Figure  2-5  for  quiet  and 
frontal  conditions.  For  the  quiet  data,  we  see  that  the  surface  for  most 
rf  and  envelope  levels  shows  no  correlation;  it  appears  as  a  product 
of  an  envelope  type  density  and  an  rf  density.  At  large  envelope  values, 
however,  the  conditional  rf  density  (a  line  of  constant  envelope  value) 
begins  to  shift  to  the  right,  indicating  larger  rf  values  are 
more  probable  at  increasing  values  of  the  pilot  channel  envelope. 

The  frontal  surface  shows  much  the  same  behavior,  except  that  the 
dependency  of  the  rf  density  on  the  conditioning  envelope  value  extends 
over  nearly  all  of  the  observed  range  of  these  variables.  In  Fig¬ 
ure  2-6,  the  joint  probability  surface  for  tropical  conditions  is  shown, 
and  exhibits  a  dependency  behavior  intermediate  between  the  quiet  and 
frontal  surfaces.  Figure  2-6  also  shows  the  one  data  recording  made 
of  joint  envelope  dependency  between  two  frequency  channels.  The  move¬ 
ment  of  the  peak  of  the  envelope  density,  parallel  to  the  conuitioning 
envelope  axis  and  then  to  the  right  for  larger  values,  provides  a  clearer 
indication  of  the  nature  of  the  statistical  interrelation  between  noise 
waveforms  in  disjoint  channels. 
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A  more  intuitive  measure  of  this  interrelationship,  that  can  be 
displayed  in  two  dimensions,  is  the  standard  deviation  of  the  rf  noise 
waveform,  conditioned  on  the  value  of  the  envelope  in  the  pilot  channel. 
For  our  discrete  sample  record  this  is  given  as 
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Plots  of  this  joint  channel  statistic  for  various  noise  conditions  and  band¬ 
width  are  shown  in  Figure  2-7.  These  plots  were  computed  directly  from 
the  raw  data,  not  smoothed  estimates  used  in  the  joint  density  surfaces. 
Note  that  the  relative  scale  factors  used  for  different  data  records, 
removed  in  Figure  2-1,  are  present  here.  The  horizontal  portion  of  each 
of  the  plots  in  Figure  2-7  corresponds  to  the  region  of  the  joint  density 
surface  where  there  is  no  statistical  dependence  between  channels.  Phys¬ 
ically,  this  portion  is  due  to  the  background  atmospheric  noise  component 
where  one  cannot  identify  the  effect  of  any  single  discharge.  The  rising 
portion  of  the  plot  corresponds  to  those  components  of  the  noise  wave¬ 
form  that  are  clearly  influenced  by  a  dominant  lightning  discharge,  and 
a  linear  dependence  is  seen  in  this  joint  channel  statistic.  It  is  also 
important  to  note  that  this  linear  dependence  holds  for  all  noise  condi¬ 
tions,  from  quiet,  where  only  10%  of  the  samples  fall  in  the  correlated 
region,  to  frontal,  where  90%  fall  in  this  region/  The  nature  of  the 
joint  channel  noise  dependence  will  be  considered  in  greater  detail  in 
Chapter  3  in  connection  with  the  mathematical  noise  model. 


'  The  method  used  to  determine  these  percentages  is  given  in  Chapter  3. 
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2.  5  Time  Structure  of  Atmospheric  Noise 
2.  5. 1  Power  Spectral  Density 

g 

A  large  number  of  experimental  observations  show  that  broadband 
atmospheric  noise  has  an  approximate  l/f  dependence  from  10  kHz  to 
200  MHz.  When  such  noise  is  observed  through  a  narrow-band  filter, 
we  would  expect  that  the  noise  at  the  filter  input  may  be  considered  to 
have  a  flat  power  density  spectrum  (white  noise)  across  the  effective 
filter  bandwidth.  Thus  the  power  density  spectrum  of  the  bandlimited 
noise  should  be  proportional  to  the  filter  response  magnitude  squared, 
or  equivalently,  the  autocorrelation  of  the  bandlimited  noise  should  be 
proportional  to  the  autocorrelation  of  the  filter  impulse  response. 

To  test  this  hypothesis,  the  envelope  of  the  autocorrelation  coefficient 
of  bandlimited  atmospheric  noise  was  measured  by  sampling  the  wave- 
lorm  at  an  integral  multiple  of  the  period  of  the  nominal  center  frequency 
of  the  narrow-band  filter  and  autocorr elating  this  sample  record.  If  one 
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uses  Rice's  representation  for  a  bandlimited  waveform,  it  is  well  known 
that  the  autocorrelation  function  is  given  as 

R(t)  =  Raa(x)  cos  (u>ot)  +  Rab(r)  sin  (wqt), 

where  R, .  (t)  is  produced  by  the  power  spectral  density  component  sym- 
metric  with  respect  to  the  frequency,  wq,  and  R^t)  by  the  anti -symmetric 
component.  By  choosing  the  arbitrary  center  frequency,  u>q,  as  the  fun¬ 
damental  of  our  sampling  frequency  (nominally  band  center)  we  see  that 
the  autocorrelation  of  our  sample  record  will  be 

R(nT)  =  R  (nT), 

cla 
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which  is  the  symmetric  component  of  the  autocorrelation  envelope. 

This  component  of  the  autocorrelation  envelope  was  estimated, 
using  an  8000-point  sample  record,  for  three  different  filter  shapes. 

The  filter  was  alternately  excited  by  atmospheric  noise  and  a  known 
white  gaussian  signal  source  (General  Radio  Type  1390-B  Random  Noise 
Generator).  The  experimental  autocorrelation  envelope  is  shown  in 
Figure  2-8.  Two  sources  of  error  affected  these  measurements,  slight 
differences  in  the  sampling  frequency,  which  produce  a  slightly  differ¬ 
ent  symmetric  spectral  component,  and  the  greater  instability  of  the 
estimate  with  atmospheric  noise.  The  first  factor  produces  slightly 
different  shapes  of  the  main  lobe  at  the  origin,  while  the  latter  produces 
greater  fluctuations  in  the  "long-time"  correlation  estimates.  The 
results  shown  in  Figure  2-6  appear  to  confirm  the  general  hypothesis 
that  received  atmospheric  noise  fields  may  be  considered  white  with 
respect  to  typical  filter  bandwidths.  The  import  of  this  fact  to  signal 
processing  problems  is  that  no  unique  information  is  conveyed  by  knowl¬ 
edge  of  the  bandlimited  atmospheric  noise  correlation  function,  beyond 
that  provided  by  a  white  noise  assumption. 

2.  5.  2  Time  Structure  of  A(t) 

The  second  most  important  aspect  of  our  experimental  investigation, 
after  the  measurements  of  joint  channel  characteristics,  was  an  estimate 
of  the  time  structure  (dynamics)  required  of  the  lowpass  modulating  ran¬ 
dom  process,  A(t),  in  the  multiplicative  noise  model  y(t)  =  A(t)  •  n(t). 

The  integrated  estimates,  A(t.),  were  used  to  explore  this  behavior. 
Preliminary  tests  of  these  estimators  (see  Appendix  C)  indicated  that 
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tho  1  and  10  msec  integration  periods  provided  a  low  variance  estimate 
of  A(t.)  when  used  with  a  gaussian  signal  source  and  rf  bandwidths  of 
1  and  10  kHz,  and  clearly  exposed  a  periodic  two-state  power  level  used 
to  modulate  the  gaussian  source. 

Recordings  of  the  A(t.)  estimates  were  made  for  all  noise  conditions 
at  frequencies  of  14,  65,  and  83  kHz.  The  100  kHz  channel  could  not  be 
used  since  these  A(t.)  records,  used  mainly  for  correlation  analysis, 
were  made  in  continuous  bursts  of  0.  5  sec  and  longer  and  such  lengths 
would  have  been  contaminated  by  the  pulsed  signal  groups  of  the  operating 
Loran-C  navigation  system.  The  A(t.)  samples  were  taken  simultane¬ 
ously  in  two  frequency  channels,  with  the  1  msec  integration  periods 
used  to  study  the  short-time  structure  of  A(t),  and  the  10  msec  estimates 
the  long-time  structure.  Physically,  this  corresponded  to  emphasizing 
the  structure  of  single  bursts  (average  time  of  1  msec)  and  multiple  dis¬ 
charges  (average  time  300-400  msec). 

Prior  to  a  discussion  of  the  correlation  analysis  of  these  A(t.) 
records,  it  is  of  interest  to  inspect  an  actual  sample  record  as  shown 
in  Figure  2-9.  This  is  a  plot  of  each  sample  point  (the  digitized  inte¬ 
grator  output  at  1  and  10  msec)  connected  by  a  straight  line,  with  the 
appropriate  time  scale  indicated.  We  see  from  the  upper  1  msec  inte¬ 
gration  time  the  tendency  of  the  noise  fluctuations  to  cluster  in  bursts, 
and  also  the  greater  fluctuation  of  the  65  kHz  estimates  in  comparison 
with  the  14  kHz  estimates.  This  is  due  to  the  intense  leader  structure 
feund  at  these  frequencies  and  above.  The  10  msec  sample  record 
makes  this  behavior  more  apparent,  where  there  is  a  somewhat  lower 
range  of  excursions  in  the  A(t.)  estimates  in  the  65  kHz  trace,  but  a 


Figure  2-9.  Sample  Record  of  A(t.)  Estimates 
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larger  fluctuation  in  the  continuum  component.  This  undulating  behavior 
is  suggestive  of  a  lowpass  modulating  process,  with  the  addition  of  a 
larger  white  or  uncorrelated  component. 

The  autocorrelation  of  the  estimates  A(t.)  with  1  msec  integration 
times  showed  the  same  basic  structure  under  all  noise  conditions  as 
shown  in  Figure  2-10.  The  A(t^)  time  series  is  dominated  by  an  uncor¬ 
related  component  producing  the  large  discrete  impulse  at  the  origin, 
followed  by  a  small  exponential  type  short  time  decay  component.  These 
records  were  all  made  with  a  1  kHz  bandwidth;  however,  previous 
records  using  a  10  kHz  bandwidth  showed  the  same  behavior.  We  attrib¬ 
ute  this  short-time  decay  to  the  stochastic  length  of  single  discharges, 
including  the  leader  stroke.  The  difference  between  the  value  of  the 
autocorrelation  coefficient  at  3-4  msec  and  that  at  infinity  (computed 
as  the  square  of  the  average  value  of  A(t.)>  is  caused  by  the  presence 
of  long  multiple  discharges,  especially  during  the  tropical  and  frontal 
conditions. 

The  autocorrelation  coefficient  of  the  10  msec  estimates,  A(t), 
show  the  exponential  type  decay  due  to  this  long-time  feature  as  seen 
in  Figure  2-10.  This  figure  shows  the  discrete  impulse  at  the  origin, 
caused  by  the  same  impulse  measured  with  the  1  msec  estimates  plus  the 
short-time  correlation  component  observed  there,  followed  by  the  long¬ 
time  decay  which  has  a  time  constant  in  the  300-500  msec  category.  To 
study  this  decay  component  several  recordings  were  made  with  the  start 
of  a  sample  burst  triggered  by  the  10  msec  estimate,  A(t),  exceeding 
a  preset  threshold.  The  10  msec  estimate  would  only  exceed  this 
threshold  when  a  strong  multiple  discharge  had  occurred,  sufficient  to 


Figure  2-10.  Short-Time  and  Long-Time  Autocorrelation  of  A(t.) 
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cause  A(t.)  to  rise  well  above  the  level  normally  caused  by  the  background 
noise  plus  the  average  single  discharge  rate.  The  autocorrelation  of  this 
10  msec  conditional  sample  record  for  tropical  conditions  is  shown  in  Fig¬ 
ure  2-11  where  the  exponential  decay  is  strongly  in  evidence.  We  also  note 
the  dramatic  difference  in  magnitude  of  the  decay  component  between  14kHz 
and  65  kHz  due  to  the  difference  in  the  lightning  burst  structure  at  these 
frequencies.  Due  to  the  effect  of  the  conditioning  circuit  in  these  observa¬ 
tions  we  should  not  draw  general  conclusions  that  this  large  a  difference 
is  always  present,  as  can  be  seen  from  the  sample  record  in  Figure  2-9. 

During  the  quiet  noise  period  when  there  was  only  scattered,  rela¬ 
tively  low  energy,  discharge  activity,  there  was  little  significant  multiple 
discharge  phenomena.  This  is  shown  in  the  autocorrelation  of  the  10  msec 
sample  record  from  this  period  (Figure  2-12)  in  which  there  is  little  dis¬ 
cernible  long-time  decay.  The  autocorrelation  of  the  10  msec  sample 
record  for  quiet-night  conditions  (Figure  2-12)  does  show  the  long-time 
decay  factor.  This  was  caused  by  the  north-moving  warm  air  soon  to 
cause  locally  tropical  noise  conditions,  and  the  fact  that  propagation  dis¬ 
tances  are  greatly  extended  at  night  by  ionospheric  reflections,  thus 
extending  the  region  of  geographic  noise  sensitivity  towards  the  southern 
warm  air  mass. 

The  various  cross -frequency-channel  correlations  shown  in  Fig¬ 
ures  2-10  through  2-12,  which  reach  values  of  0.95  in  some  cases, 
reinforce  the  modeling  concept  of  A(t)  being  frequency  independent.  We 
note  that  cross-correlation  is  generally  better  for  the  65  kHz-83  kHz 
pair  than  for  the  65  kHz- 14  kHz  pair.  This  is  to  be  expected,  due  both 
to  the  increased  frequency  separation  and  to  the  changing  nature  of  the 
dominant  discharge  mechanism  at  the  14  kHz  frequency. 


Figure  2-11 
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The  principal  conclusion  that  we  draw  from  the  correlation  analysis 
of  the  A(L)  sample  records  is  that  there  is  not  a  significant  time  struc¬ 
ture  in  the  time-varying  noise  power  level.  By  significant,  we  mean  a 
structure  that  can  be  exploited  in  receiver  design  to  estimate  future 
values  of  the  power  level  from  present  observations.  We  shall  find  in 
Chapter  4  that  the  noise  t'me  structure  does  have  a  demonstrable  effect 
on  receiver  performance  and  that  this  can  be  largely  compensated  for 
by  use  ol  the  joint  channel  information  to  estimate  directly  both  the 
correlated  and  uncorrelated  component  of  A(t),  without  consideration 
of  a  model  for  A(t)  dynamics.  Our  experimental  A(L)  measurements 
will  be  used  in  the  noise  model  development. 

2. 6  Summary 

In  this  chapter  we  have  presented  a  number  of  statistical  measure¬ 
ments  of  the  characteristics  of  bandlimited  VLF  and  LF  atmospheric 
noise  waveforms.  These  characteristics  were  measured  and  interpreted 
from  the  point  of  view  of  a  postulated  multiplicative  noise  model,  y(t)  = 
A(t)  *  n(t).  The  important  results  may  be  summarized  as  follows. 

a)  Atmospheric  noise  waveforms,  especially  at  LF,  exhibit  a  wide 
variation  in  characteristics  dependent  upon  local  weather  condi¬ 
tions.  Since  the  two  physical  factors,  propagation  and  noise 
mechanisms,  become  more  important  at  frequencies  above  LF, 
we  expect  this  variation  is  similarly  more  pronounced  at  higher 
frequencies. 

b)  Measurement  of  first  order  joint  statistics  between  a  bandlimited 
rf  noise  waveform  in  frequency  channel  one  and  an  envelope 
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from  frequency  channel  two  show  a  significant  statistical  depen¬ 
dence  above  a  threshold  value.  The  model  concept  that  A(t)  is 
independent  of  frequency  attempts  to  describe  these  observations. 

c)  The  received,  broadband,  atmospheric  noise  field  appears  as 
an  uncorrelated  random  process  when  observed  via  narrow-band 
filters.  The  average  power  spectral  density  of  the  observed 
narrow-band  noise  is  proportional  to  the  square  of  the  magnitude 
of  the  filter's  frequency  response. 

d)  The  hypothetical  modulating  process,  A(t),  should  be  basically 
uncorrelated  beyond  1  msec,  with  only  a  small  short-time  and 
long-time  dynamic  damping.  This  observation  is  only  valid  for 
observation  bandwidths  of  1  kHz  or  larger. 

In  the  next  chapter  we  shall  develop  a  noise  model  based  on  the  mul¬ 
tiplicative  concept  which  fits  the  observed  data  over  the  wide  range  of 
observed  noise/ weather  conditions. 


Chapter  3 

AN  ATMOSPHERIC  NOISE  MODEL 
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In  Chapter  1  of  this  report  we  gave  the  motivational  background  for 
a  multiplicative  model  for  bandlimited  atmospheric  radio  noise.  Chap¬ 
ter  2  presented  a  number  of  experimental  observations  of  characteristics 
such  a  model  must  possess.  In  this  chapter,  we  shall  develop  the  spe¬ 
cific  mathematical  form  of  a  multiplicative  model.  This  model  will  be 
compared  with  the  observed  data  and  it  will  be  shown  that  the  necessary 
model  parameters  can  be  chosen  in  a  logical  manner.  Principal  weight 
will  be  given  to  the  modeling  of  the  first  order  joint  channel  statistics 
since  the  observations  of  Chapter  2  show  that  these  are  a  distinguishing 
characteristic  of  atmospheric  noise.  A  method  of  generating  an  approxi¬ 
mation  to  the  time  structure  of  A(t)  will  then  be  given  and  it  will  be  shown 
that  using  the  parameters  developed  for  the  first  order  model,  acceptable 
simulation  results  of  the  A(t.)  estimates  are  obtained.  Finally,  a  can¬ 
onic  atmospheric  noise  generator,  suitable  for  computer  implemen¬ 
tation  using  uniform  and  gauss ian  random  number  generators  will  be 
given.  This  is  accompanied  by  a  table  specifying  values  of  the  model 
parameters  which  describe  the  range  of  experimental  data  reported  in 
Chapter  2. 

3. 1  Hall's  Noise  Model 
1 8 

Hall  suggested  the  multiplicative  noise  model  y(t)  =  A(t)  n(t)  as  a 
description  of  a  bandlimited  atmospheric  noise  waveform  in  a  single¬ 
frequency  channel.  Hall  proposed  that  the  lowpass  random  process  was 
generated  by  the  inverse  of  a  chi  process,  that  is. 
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A(t)  = 


m 

2 

i=l 


(3.1) 


where  the  b.(t)  are  statistically  independent  gaussian  (N(0,o-))  lowpass 
processes.  The  first  order  pdf  of  A(t),  the  inverse  of  a  x(m,  <r)  ran¬ 
dom  variable,  is 


exp 


-m 


h*2A2 


-oo  <  A  <  +oo.  (3.  2) 


From  (3.  2)  and  the  assumption  that  n(t)  is  gaussian  (N(0,crnjj)  we  can 
compute  the  pdf  of  y(t)  as 
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(3.3) 


where  we  have  set  the  cr  parameter  of  the  chi  process  equal  to  1  for 
convenience  since  it  only  appears  as  a  divisor  of  cr-j.  Recalling  the 
measured  atmospheric  noise  pdf's  given  in  Figure  2-1,  we  see  that 
expression  (3.  3)  has  the  correct  basic  shape  to  fit  these;  a  slope  of 
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2  2 

zero  as  y  —  0,  a  single  "break-point”  on  a  log-log  plot  at  y  >  mo-^  and 
a  hyperbolic  decay  beyond  that  point,  appearing  as  straight  line  of 
slope  -(m+1)  on  a  log  vs  log  plot. 

Expression  (3.  3)  is  a  valid  pdf  for  any  m  >  0,  although  Hall  generally 
restricted  m  to  be  integer.  A  best-fit  selection  of  m  can  be  made 
from  the  large  amplitude  region  of  the  data  pd .  Hall  also  suggested  that 
the  dynamics  of  A(t)  could  be  provided  by  specification  of  a  ratirina]  a 
power  spectral  density  for  b^(t).  Thus  we  see  that  m  has  a  second  im¬ 
portant  role  in  the  degree  of  freedom  allowed  to  effect  the  dynamics  of 
A(t). 

Hall  provides  a  similar  closed  form  expression  for  the  atmospheric 
noise  envelope,  from  which  he  v.  omputes  the  exceedance  probability  (one 
minus  the  cumulative  probability)  of  the  envelope  to  compare  with  pub- 

7 

lished  data,  principally  that  given  by  Watt.  He  concludes  that  his  model 
provides  reasonable  fits  to  these  data,  with  a  value  of  m  =  2  fitting  most 
data  reported  for  moderate  mid-latitude  noise  conditions  at  VLF.  Hall 
also  computes  a  level  crossing  statistic.  However,  since  we  have  taken 
a  different  experimental  approach  to  considering  time  structure,  we  shall 
not  deal  further  with  the  latter  comparisons. 

3.  2  An  Extension  of  Hall's  Model 

If  we  compare  the  pdf  of  (3.  3)  in  greater  detail  with  the  observed 
probability  densities  of  Figure  2-1,  we  see  that  (3.3)  does  not  have  suf¬ 
ficient  degrees  of  freedom  to  account  for  the  second  break  point  seen  in 
the  tail  of  the  quiet  noise  pdf.  The  occurrence  of  this  breakpoint  is  indi¬ 
cative  of  the  noise  waveform  becoming  closer  to  a  gaussian  process. 
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and  by  implication  our  model  should  have  this  capability  of  being  reduced 
to  a  gaussian  process.  A  similar  difficulty  is  encountered  when  we  con¬ 
sider  the  joint  channel  first  order  data  (Figures  2-5  and  2-6).  The  gen¬ 
erating  mechanism  of  (3. 1)  does  not  appear  to  provide  for  the  uncorre¬ 
lated  portion  of  the  conditional  rf  noise  standard  deviation,  if  we  make 
our  assumption  that  A(t)  is  deterministically  identical  for  both  channels. 

A  somewhat  less  important  difficulty  is  encountered  if  we  attempt 
to  associate  a  time  structure  for  the  b^(t)  as  that  provided  By*a  linear  **' 
lowpass  filter  excited  by  white  gaussian  noise.  As  was  seen  m  Fig¬ 
ures  2-10  through  2-12  the  modulating  process,  A(t),  must  be  domi¬ 
nated  by  an  uncorrelated  component,  but  have  two  second  order  time 
correlations.  However,  as  will  be  seen,  m  values  of  2  or  1  are  re¬ 
quired  to  match  our  observations  and  these  do  not  provide  sufficient 
freedom  in  the  A(t)  generator  mechanism  of  (3. 1)  to  meet  these  time 
structure  observations  of  Chapter  2. 

The  basic  difficulty  noted  in  the  first  order  behavior  of  (3. 1)  is 
that  the  bandlimited  atmospheric  noise  waveform  has  a  distinct  period 
when  it  is  in  a  nearly  gaussian  state,  the  so-called  background  compo¬ 
nent  we  have  referred  to.  Our  first  attempts  to  deal  with  this  in  terms 
of  A(t)  consisted  of  adding  a  constant  bias  term  to  the  dynamic  A(t)  gen¬ 
erator  of  (3. 1).  This  has  the  effect  of  suppressing  the  dynamic  compo- 

19 

nent  until  it  exceeds  a  threshold,  and  this  idea  leads  naturally  to  Kapp's 
idea  of  a  switched  process  model.  We  have  adapted  this  to  Hall’s  model 


as  our  proposed  model,  given  as 
y(t)  =  nj(t)  +  x(t)  a(t)  n2(t), 


(3.4) 
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where  y(t)  is  a  bandlimited  atmospheric  noise  waveform.  The  process 
a(t)  ng(t)  is  identical  to  that  given  by  Hall.  We  shall  refer  to  it  as  the 
"Hall  component"  and  use  for  the  associated  parameter  of  ngCt). 
x(t)  is  a  two-state  process,  assuming  values  one  or  zero  with  probability 
Px  and  qx,  respectively.  The  process  n^(t)  is  a  narrow-band  gaussian 
process  (first  order  pdf  N(0,  o^))  and  provides  the  background  gaussian 
component  required  of  our  model.  The  components  n^t)  and  ng(t)  are 
» »-»etatistically  independent,  with  power  spectral  densities  proportional  to 
the  bandlimiting  filter's  magnitude  response  squared.  The  function  of 
x(t)  is  to  "switch-on"  the  non-gaussian  or  Hall  noise  component,  mod¬ 
eling  those  time  periods  when  the  atmospheric  noise  waveform  is  in  a 
distinctly  non-gaussian  state.  In  addition  to  providing  several  more 
degrees  of  freedom  for  first  order  statistics,  we  also  have  available  the 
time  structure  of  x(t),  coupled  with  that  of  a(t),  to  model  the  observations 
of  A(t)'s  time  structure. 

The  first  order  pdf  of  y(t),  defined  by  (3.4),  is  given  by 
Py(y)  =  N(0,<rnG)  ®  [q6(r1)  +  pxpyH(T1)] 

=  qN(0,  o-nG)  +  PxN(0,  <rnG)  ®  Py^n),  (3.  5) 

where  ®  denotes  convolution,  6(*)  the  Dirac  delta  function,  and  Pyjj(T)) 
is  the  Hall  pdf,  (3.3).  The  characteristic  function  (cf)  of  p  (y)  can  be 

y 

expressed  in  closed  form,  but  we  have  not  been  able  to  obtain  its  inverse 
transform.  We  note  that  the  background  process,  n^(t),  is  quite  small 
with  respect  to  the  large  noise  amplitudes  provided  by  the  Hall  compo¬ 
nent,  hence  the  gaussian  pdf  appears  as  an  impulse  scanning  function  at 


these  large  values  of  the  argument  y.  Thus  our  model  preserves  the 

1 

hyperbolic  tail  (Lim  p  (y)  cc  l/y  )  of  the  Hall  component  to  match 
y-00  y 

the  observed  large  amplitude  behavior  of  atmospheric  noise.  In  most 
practical  cases,  the  resulting  p  (y)  can  be  approximated  as  a  sum  of  a 

J 

gaussian  and  Hall  pdf;  however,  we  have  used  numerical  techniques  to 
evaluate  p  (y)  in  section  3.6.2  below. 

v 

3.  3  Noise  E nvelope 

To  develop  an  approximation  to  the  envelope  pdf  for  a  bandlimited 

noise  waveform  given  by  (3.  4),  we  begin  with  the  definition  of  the  joint 
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envelope  and  phase  pdf  as 


2  2  v  2 

p  Av,  <)>)  =  vp  (v  cos  <J>,  v  sin  <)>),  v  =  y  +  y  , 

y*y 

()>  =  tan  y/y,  y  -  Hilbert  Transform  (y). 


(3.6) 


In  our  case 


y  =  n,  +  xan9. 


V  y  V 

y  =  nj  +  xan2. 


y 

where  for  2  gaussian  narrow-band,  n^  2  are  also  gaussian  narrow- 
band  and  statistically  independent.  The  quantity  xa  need  not  be  hatted 
if  we  assume  that  xa  has  no  frequency  overlap  with  n2,  that  is,  it  is 
slowly  varying  with  respect  to  the  band -center  frequency.  Thus  the 
envelope  pdf  is  formally  given  as 


-62- 


p  .  (v,  4>)  =  v  \  \  p  (v  cos  <J>,  v  sin  <i> ,  a,  x)  dadx 
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The  last  integral  with  respect  to  "a"  is  similar  in  form  to  that  leading 
to  (3,  3)  and  is  recognized  as  Hall's  expression  for  the  envelope  pdf  of 
his  noise  model.  Expression  (3.  7)  suggests  a  lorm  of  the  first  order 
envelope  pdf  (with  phase  random  0-2tt)  similar  to  that  given  for  the  band- 
limited  noise  by  (3.  5),  if  we  identify  the  second  term  as  a  convolution  of 
a  Rayleigh  pdf  with  the  Hall  envelope  pdf.  The  second  term  of  (3.7)  does 
not  reduce  to  this  convolution;  however,  we  shall  use  this  form  as  an 
approximation  to  the  noise  envelope.  This  approximation  then  yields  an 
envelope  and  associated  first  order  pdf; 

v(t)  =  Vj(t)  +  x(t)  a(t)  VgU),  (3.8) 

where  v^(t)  is  Rayleigh  (R(o-q))  and  Vg(t)  is  Rayleigh  (Rfcr^)),  the  pdf  of 


We  note  that  this  approximate  model  and  resulting  pdf  for  the  envelope 

provides  the  correct  small  amplitude  behavior  (Lim  p  (v)  oc  v)  and  a  large 

/V-0  V  1  \ 

amplitude  behavior  following  a  hyperbolic  tail  (Lim  p  (v)  oc  — — -  1. 

v—oo  v  vm 

3. 4  Joint  Channel  Model 

We  can  now  combine  the  model  for  a  bandlimited  atmospheric  noise 
waveform,  observed  in  frequency  channel  fl,  the  approximate  envelope 
model  for  frequency  channel  f2,  and  our  assumption  that  a(t)  is  identical 
in  both  models  to  determine  joint  channel  characteristics.  When  x(t)  =  0, 
or  the  Hall  noise  component  is  very  small  (a(t)  small)  the  background 
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gaussian  noise  component  dominates  both  t  .velope  and  rf  noise  densities 


and  the  joint  channel  density  is  simply  the  product  of  the  two  since  they 
are  statistically  independent.  When  x(t)  =  1  and  the  Hall  component  is 
dominant,  we  can  proceed  exactly  as  in  developing  the  envelope  and  rf 
Hall  noise  densities  ((3.  3)  and  (3. 9))  to  develop  the  joint  channel  pdf. 
Thus 


Py.  v** V)  =  PyH,  vH(y-  v>  =  J  Py,  v  |  V'  4  PaW  d,> 

9- 

*  J  Py|a(y-,')pv|a(v',l)pa<,')dT' 


'I 


00 


0  T>%H'ra' 

2(m/2)m/2 


exp 


4 

-y 

V 

2 

-v 

,2  2 
2ffnH^ 

2  2  exp 
*1  ^vH 

i.2  2 

2%H1  J 

■(f)' 


m+1  expl  2r|2  Jdrj 


I n  m+1 

(f ) Tisr-) 2  2  v 


(JL)2+(^_/  +  m 

V<rnH;  V<rvHy 


m+3  » 
2 


-oo  ^  y  <  +oo,  0  ^  v  <  +oo.  (3.  1 1) 


In  a  similar  manner,  the  joint  pdf  for  the  Hall  component  of  the  envelope 
in  two  disjoint  frequency  channels  can  be  computed  as 
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>  v  <viJv2)  = 
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(3.  12) 
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From  the  joint  pdf's  and  the  single -channel  pdf's  we  can  determine  the 
conditional  pdf  of  either  variable.  For  example,  the  conditional  pdf  of 
the  rf  noise,  y,  given  a  pilot  channel  envelope  value  v  ,  is 


P  (y,  v  ) 
y,  v  J *  o 
o 
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p  (v  ) 
*v  o 
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oc 


,  m+3 
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2  +  ' 

TnH  °vH 


+  m 


(3.13) 


Comparing  (3. 13)  to  the  unconditional  pdf  of  y,  (3.3)  and  the  observed 
data  in  Figures  2-5  and  2-6  we  see  that  it  has  the  correct  form  for 
small  values  and  large  values  of  y,  and  that  the  conditional  pdf  tail 
decreases  by  a  power  of  two  faster  than  the  unconditional  pdf.  This 
behavior  will  be  considered  in  more  detail  in  section  3.6.  3  below. 


3.  5  First  Order  Model  Statistics 
3.  5. 1  Unconditional  Moments 

Inspection  of  the  limiting  large  amplitude  behavior  of  the  Hall  com¬ 
ponent  of  the  noise  pdf's  ((3.  3  and  (3.  9)),  given  by 
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(3. 14b) 
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shows  that  both  have  unbounded  first  and  second  moments  for  m  =  1  and 
unbounded  second  moments  for  m  =  2.  These  are  the  values  of  m  re¬ 
quired  to  fit  all  of  our  observed  data  and  we  are  thus  faced  with  a  situa¬ 
tion  of  attempting  to  model  a  physical  phenomenon  with  a  model  that  has 
an  infinite  power  level.  Hall  suggested  this  difficulty  be  overcome  by 
truncating  the  range  of  the  inverse  chi  process,  a(t).  We  shall  simply 
truncate  the  complete  noise  density  at  a  value  of  y  or  v  which 

IUcLX  I XI  cl  X 

corresponds  to  the  upper  limit  of  our  observed  noise  density.  This 
truncation  was  the  most  convenient  in  our  numerical  comparisons  of  the 
model  to  the  observed  noise  characteristics.  In  the  case  of  realizing 
the  canonic  atmospheric  noise  generator  (to  be  given  in  section  3.8)  in 
a  Monte  Carlo  computer  simulation  it  is  probably  more  convenient  to 
truncate  a(t).  Since  ng(t)  is  gaussian,  the  resulting  distribution  will 
decay  rapidly  beyond  such  an  a(t)  truncation  point  and  hence  it  is  prob¬ 
ably  not  of  great  practical  significance  which  method  is  employed  since, 
as  we  shall  find,  practical  signal  processing  structures  are  quite  insen¬ 
sitive  to  the  exact  probability  of  the  large  amplitude  noise  components. 

Distribution  moments,  the  second  moment  or  power  level  in  par¬ 
ticular,  are  well  ingrained  in  engineering  thinking  due  to  their  role  in 
linear  systems  and  gaussian  noise.  However,  for  non-gaussian  noise, 
they  are  only  two  parameters  of  a  theoretically  infinite  set  required  to 
completely  characterize  an  arbitrary  distribution.  In  the  case  of  non- 
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gaussian  atmospheric  noise,  they  are  only  moderately  informative,  as 
when  the  ratio  of  rms  to  mean  value  of  the  envelope,  the  so-called 
parameter,1^’  is  used  to  characterize  the  degree  of  non-gaussianity 

of  the  pdf.  We  shall  find  in  Chapter  4  that  a  simple  nonlinear  system 
will  perform  upwards  of  17  dB  better  in  frontal  type  of  noise  as  com¬ 
pared  with  gaussian  noise  of  equal  power  level.  Improved  signal  pro¬ 
cessing  structures  of  this  type  are  sensitive  to  the  nature  of  the  noise 
pdf  in  the  low-to-middle  amplitude  region  and  almost  insensitive  to  the 
large  amplitude  portion,  while  the  latter  in  many  cases  can  dominate  the 
second  moment  or  power  level.  For  these  reasons,  as  well  as  the  prac¬ 
tical  one  that  measurement  of  the  noise  power  level  with  analog  instru¬ 
ments  is  a  very  difficult  proposition  due  to  the  dynamic  range,  we  have 
considered  first  and  second  moments  of  the  noise  model  only  incidentally, 
after  the  model  parameters  have  been  selected  by  other  criteria. 

Values  of  the  u~st  and  second  moments  of  the  rf  and  envelope  noise 
model  pdf's  are  given  below  for  the  cases  of  interest. 

a)  General 

|y|  =  I  nj+xangl  «  |nj|  +  xty^ 
v  -  Vj  +  xavg  =  v’j  +  x~H 

|  y|  2  «  n2  +  2x|yJ|nJ+ x2  y2 

v2  =  v2  +  2vl  vHx  +  x2  v2  . 

Using  Gaussian/Rayleigh  and  Hall  component  pd.'^  with  upper  integration 


(3.22) 


3.  5.  2  Conditional  Moments 
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Interchanging  the  order  of  integration  and  substitution  of  a  gaussian  pdf 

for  py*i|a<T1’^  yields 


1 
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Substitution  of  a  Rayleigh  pdf  for  PvH|a(v0*  D  and  the  pdf  ^3‘ yields 
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(3.23) 


These  two  results  are  conditioned  on  both  the  value  of  the  pilot  chan 
nel  envelope,  v  ,  and  the  state  of  x(t),  whereas  the  data,  such  as  Fig¬ 
ure  2-7,  is  conditioned  only  on  vQ.  Thus,  these  results  are  valid  for 
small  and  large  amplitudes  of  the  envelope,  where  one  may  infer,  in  a 
Bayesian  sense,  the  state  of  x(t)  being  0  or  1,  respectively.  These 
results  are  sufficient  to  describe  all  but  a  small  transition  region  of  the 
conditional  rf  noise  variance. 

Using  the  approximation  for  the  noise  envelope  (3.8)  and  the 
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Rayleigh  and  Hall  component  pdf's  we  obtain  similar  results  for  the 
conditional  second  moment  of  the  joint  envelope -pilot  channel  envelope; 

a)  with  x(t)  =  0 


Jv  Ar1|v2,x=0(,l'v2)d,'=  2°Ig- 


(3.24) 


b)  with  x(t)  =  1, 

1  'l2pv1H|v2H,x-0(,'-v2)  d’>  = 


2  tr 


v,H  r 


v2H 


v2  ■"  m,rv2H. 


(3.25) 


We  shall  also  compare  the  ratio  of  mean  to  variance  of  the  condi¬ 
tional  | rf |  channel  pd  for  the  model  and  the  data,  as  a  measure  of  the 
non-gaussian  character  of  the  conditional  rf  distribution.  This  is  com¬ 
puted  from  the  joint  y-v  pdf  (3. 12),  with  m  =  1,  as  follows.  The  use  of 
the  upper  limit  of  integration  is  not  required  here  because  the  moments 
are  unbounded,  rather  it  includes  the  effect  of  the  maximum  range  of 
the  observations  which  affects  the  ratio  of  interest  at  large  values  of 
the  envelope. 
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3. 6  Comparison  of  Joint  Channel  Model  with  Data 
3.6. 1  General  Parametric  Dependence 

We  have  used  the  unconditional  probability  densities  of  the  rf  noise 
from  frequency  channel  1  and  the  envelope  from  frequency  channel  2, 
plus  the  conditional  standard  deviation  of  the  rf  noise,  as  the  data  base 
to  compare  our  first  order  model  with  observed  data.  The  joint  chan¬ 
nel  model,  with  parameters,  is  given  by 


nG  nH 

yfi  =  n|(t)  +  x(t)  a(t)  n2(t) 


Px  m 


vf2  =  vj(t)  +  x(t)  a(t)  v2(t). 

I  I 

°"vG  °vH 


(3.28) 
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The  parameter  m  determines  the  large  amplitude  behavior  of  the  rf  and 
envelope  pdf's  (3. 14a,  b)  and  must  satisfy  m  >  0.  We  have  restricted  m 
to  be  an  integer  which  provides  the  simple  generating  mechanism  for  the 
canonic  atmospheric  noise  generator  to  be  given  in  section  3.  8.  The 
parameter  m  is  thus  estimated  from  the  slope  of  the  observed  density 
tail  on  a  log -log  plot  (see  section  3. 1). 

We  found  in  section  3.  5.2  that  the  conditional  rf  noise  variance 

is  given  by  the  variance  of  the  background  gaussian  component,  nj  — * 

2 

°"nG*  f°r  small  values  of  the  conditioning  envelope.  This  variance  is 
then  the  square  of  the  horizontal  portion  of  the  conditional  standard  devi¬ 
ation  plot  as  seen  in  Figure  2.7.  The  parameter  p  ,  the  probability 
that  the  two  state  process  x(t)  is  in  state  1  (x(t)  =  1)  can  be  estimated 
from  a  simultaneous  plot  of  the  conditional  rf  noise  standard  deviation 
versus  the  conditioning  envelope,  and  the  cumulative  probability  of  the 
envelope.  The  point  where  the  conditional  standard  deviation  begins  to 
exhibit  a  dependence  upon  the  conditioning  envelope  is  the  envelope  value 
where  x(t)  must  be  in  state  1  to  introduce  the  dependence  via  a(t).  The 
intersection  of  this  envelope  value  with  the  cumulative  envelope  proba¬ 
bility  is  then  related  to  ^  as 


P  =  1 
x 


Pv  (V*v  ). 
o 


The  remaining  three  parameters,  orvQ»trnH'crvfj'  must  then  be 
chosen  to  simultaneously  provide  a  best  fit  to  the  two  unconditional 
probability  densities  and  the  predicted  relationship  between  the  condi¬ 
tional  standard  deviation  and  the  conditioning  envelope  given  by  (3.  23). 


3.6.2  Comparison  Results 
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The  general  algorithm  used  to  compare  the  model  with  the  observed 
joint  channel  data  was  the  following: 

a)  Estimate  the  parameters  m,  <rnQ*  and  p^  from  the  information 
given  above. 

b)  Utilize  the  computer  to  numerically  evaluate  the  two  pdf  convo¬ 
lutions  (given  by  (3.  5)  and  (3. 10))  as  a  function  of  these  three 
parameters  plus  the  three  remaining  parameters  to  achieve  the 
best  simultaneous  fit  to  the  three  data  plots.  This  was  an  itera¬ 
tive  procedure  in  which  we  subjectively  evaluated  the  data  fit, 
giving  principal  weight  to  the  fit  of  the  probability  densities 
through  the  middle  region  where  the  variance  of  the  experimen¬ 
tal  estimate  of  the  density  was  minimum  (see  Appendix  B). 

The  results  of  our  comparison  and  the  required  parameters  are  shown 
in  Figures  3-1  through  3-4  for  quiet  through  frontal  noise  conditions. 

We  see  in  Figure  3-1  in  the  upper  right  plot  the  selection  of  the  value 
of  p  at  the  break  point  in  the  conditional  standard  deviation  curve  and 

A 

the  corresponding  observation  that  the  noise  is  in  a  non-gaussian  state 
for  only  10%  of  the  observed  sample  record.  In  the  lower  right  plot  of 
the  unconditional  envelope  we  note  a  major  deviation  of  the  observed 
data  from  the  predicted  small  signal  behavior.  This  is  not  a  charac¬ 
teristic  of  the  noise  envelope  but  is  caused  by  an  instrumentation  error, 
a  dc  offset  in  the  lowpass  filter  used  to  recover  the  envelope  from  the 


The  variables  RF  and  E  are  data  variables  corresponding  to  model 
variables  y  and  v,  respectively.  Also,  p-p  =  px  and  a  jj  =  ar^. 
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wideband  linear  rectifier.  Since  this  error  does  not  appear  to  seriously 
affect  the  large  amplitude  behavior  of  the  envelope,  it  was  not  considered 
to  be  a  fatal  flaw  in  the  data.  Both  the  rf  noise  and  envelope  noise 
model  data  comparisons  show  that  the  model  does  have  sufficient  flexi¬ 
bility  to  reproduce  the  second  large  amplitude  breakpoint  which,  as  noted 
in  Chapter  2,  is  caused  by  the  low  probability-of-occurrence,  non-gaussian 
excursions  of  the  noise  waveform.  Thus  the  model  parameters  (specif¬ 
ically  px,  the  probability  that  x(t)  is  in  state  1,  adding  the  non-gaussian 
Hall  noise  component)  can  be  chosen  to  provide  a  smooth  transition  of 
the  model  to  a  purely  gaussian  noise  waveform.  On  the  conditional  stand¬ 
ard  deviation  plot,  the  line  labeled  "Model,  <r  i  "is  given  by  (3.  23), 

2  2  yfl|Vf2 

assuming  v^  >  mcr^.  The  breakpoint  predicted  by  this  relation  (i.  e. , 

where  that  assumption  does  not  hold)  occurs  approximately  at  the  actual 
data  breakpoint,  but  this  is  not  true  for  the  other  data  comparisons  so 
that  only  the  linear  portion  of  the  curve,  dependent  upon  m,  o-^  and  <rvH 
is  shown.  Finally,  the  moment  statistics  of  the  model  and  data  are  given 
in  the  title  block.  These  were  computed  from  the  relation  given  in  sec¬ 
tion  3.  5.  1,  and  we  see  that  the  model  agrees  quite  closely  with  the  data, 
as  we  would  expect  from  the  generally  good  fit  of  the  model  pdf  to  the 
data  pd. 

The  data  used  in  the  comparison  of  Figure  3-2  was  recorded  during 
a  transition  noise  condition.  The  increase  in  propagation  range  of  the 
lightning  discharge  radiation  at  night  provided  a  much  higher  geographic 
density  of  lightning  sources  which  is  reflected  in  the  px  value  of  0.  5  and 
the  disappearance  of  the  second  large  amplitude  breakpoint  in  the  pdf 
tails.  The  m  value  of  2  actually  fits  the  tails  better  here  than  in  the 
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quiet  conditions.  This  again  is  reflective  of  propagation  effects;  the 
fewer  discharges  affecting  the  quiet  data  were  local  and  hence  had  a 
tendency  to  larger  excursions  which  Implies  a  shallower  slope  or  lower 
value  of  m.  For  the  quiet-night  data,  on  the  other  hand,  the  greater  dis¬ 
tance  and  the  fact  that  more  discharges  tended  to  add  together  produced 
the  effect  of  a  steeper  slope  in  the  tail.  Again  we  note  the  generally 
excellent  match  of  the  model  with  the  composite  data  plot,  and  the  indi¬ 
vidual  channel  moments  given  in  the  title  block. 

The  remaining  two  figures.  Figures  3-3  and  3-4,  show  progres¬ 
sively  more  non-gaussian  noise  conditions  as  reflected  in  the  increasing 
value  of  px  and  the  use  of  an  m  value  of  1.  Both  of  these  data  records 
were  made  during  local  morning  or  afternoon  and  hence  tend  to  reflect 
dominance  of  individual  excursions  (x(t)  =  1)  by  strong  local  discharge 
radiation.  Multiple  discharges  are  also  prominent  during  these  condi¬ 
tions  as  seen  in  the  discussion  of  section  2.  5.  2.  We  note  that  the  mod¬ 
eling  of  the  joint  channel  dependence  continues  to  provide  good  results 
in  the  conditional  standard  deviation  plots.  The  unconditional  envelope 
model  of  Figure  3-4  is  probably  the  worst  data  match  that  we  obtained 
with  all  of  our  records;  this  same  behavior  was  found  in  both  of  the  other 
two  65  kHz  records  made  during  those  conditions.  The  dc  offset  in  the 
lowpass  filter  contributed  to  this  error,  and  it  is  also  probable  that  the 
sample  record  was  not  long  enough  (100,000  samples)  to  collect  a  repre¬ 
sentative  sample  of  the  long  multiple  discharges  that  dominate  the  frontal 
noise  conditions.  This  latter  effect  would  not  be  so  predominant  in  the 
rf  noise  record  due  to  the  randomizing  nature  of  the  phase  (or  equiv¬ 
alently  ngft))  in  this  data,  hence  the  better  match  to  the  pdf  seen  in  the 
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upper  left  plot. 

It  can  be  seen  that  the  observed  noise  conditions  span,  very  nearly, 
the  full  range  of  model  parameter  valxies.  The  frontal  conditions,  with 
p  =  0.9  and  m=  1  represent  the  largest  deviation  from  gaussianity 

.X. 

of  which  the  model  is  capable.  The  quiet  conditions,  with  p  =0.11  and 

X 

m  =  2,  approach  gaussian  noise  conditions  and  one  might  assume  that 
even  lower  values  of  px  and/or  larger  values  of  m  would  be  found  during 
temperate  latitude  winter  weather. 

A  comparison  of  the  noise  model  with  data  recorded  at  14  kHz  is 
shown  in  Figure  3-5  for  three  noise  conditions.  In  the  discussion  of 
section  3.6.  1  we  saw  the  intimate  dependence  of  the  choice  of  the  px  and 
c  G  parameters  with  the  joint  channel  dependency,  which  was  not  avail¬ 
able  for  these  VL.F  observations.  Thus  the  parameter  estimation  proce¬ 
dure  was  considerably  more  difficult,  and  more  important,  the  results 
do  not  provide  the  same  degree  of  uniqueness  in  parameter  values  as 
th-t  for  the  joint  channel  data.  For  example,  it  is  difficult  to  separate 
the  effect  of  o-^  and  p^  using  only  ihe  single  channel  data  since  these 
two  parameters  enter  in  a  product  relation  for  large  amplitudes,  as 
seen  from  the  pdf  expression  for  the  Hall  noise  component  (3.  3)  which 
is  scaled  by  p  in  the  pdf  convolution  of  (3.  5).  In  a  more  general  con- 
text,  if  we  considered  that  the  observed  probability  densities  were  three 
curves  to  be  represented  mathematically,  there  are  lour  parameters 
involved,  the  first  breakpoint,  the  following  slope,  the  second  break¬ 
point,  and  the  final  slope  (this  assumes  that  the  selected  mathematical 
expression  provides  the  correct  small  signal  behavior  and  that  pdf  nor¬ 
malization,  Py^)  dij  =  1,  provides  an  overall  scale  factor).  Since 
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our  model  provides  four  parameters,  we  cannot  infer  that  it  has  any 
unique  aspects  as  a  single  channel  model,  only  that  it  does  have  the 
flexibility  to  match  the  data  to  the  degree  shown  in  Figure  3-5. 

3.6.  3  Comparison  with  Additional  Joint  Channel  Characteristics 

To  provide  further  evidence  of  the  validity  of  the  joint  channel  model, 
the  conditional  density  of  the  rf  noise  was  measured  using  a  hardware 
conditioning  circuit.  The  rf  sample  was  taken  from  frequency  channel  1 
and  the  digitized  sample  stored  in  a  holding  register.  This  sample  was 
then  written  into  the  memory,  becoming  part  of  the  sample  record,  only 
if  the  next  envelope  sample  from  frequency  channel  2  fell  within  a  vari¬ 
able  threshold  window  determined  by  two  level  detectors.  This  technique 
allowed  longer  sample  records  at  a  given  conditioning  envelope  value 
than  were  obtained  with  the  continuous  joint  sample  records  previously 
described.  This  increased  length  provided  better  probability  density 
estimate  stability. 

If  we  assume  that  the  joint  probability  density  of  these  samples  was 
dominated  by  the  joint  Hall  component  density,  then  the  conditional  rf 
density  is  given  by 

py|  v  (y-vo! - -  '  PyH.vH^'V' 

0  pvH  vo 

Using  pdf  expressions  (3.9)  and  (3.  11)  and  transforming  to  normalized 
y'  and  v'  variables  via  the  relations  y'  =  y/ff^  and  v'  =  v/c results 
in  the  normalized  conditional  rf  noise  density  for  m  =  1 


•4- 


v' 

o 


(y',vy 


4'2+v.2+1)2' 


0  <  <  +oo,  0  <  |  y1 1  <  +oo. 

(3.29) 


Four  conditional  sample  records  were  used,  each  8000  samples  long, 
with  the  value  of  the  conditioning  envelope,  vQ,  related  between  records 
by  a  power  of  2.  The  expression  (3.  29)  was  matched  to  the  second  rec¬ 
ord  at  y’  —  0  and  at  a  large  value  on  the  tail  of  the  observed  density. 
This  two-point  fit  provided  a  value  of  vQ2n,  -1  n  2,  for  all  of  the 
records  and  the  scale  factor  for  y*  relating  the  normalized  coordinates 
to  the  observed  voltage  scales.  The  conditional  pdf  (3.  29)  was  then 
plotted  for  the  remaining  three  records  with  the  appropriate  power  of 
2  change  in  v^.  The  comparative  results  are  shown  in  Figure  3-6 
where  we  have  used  two  plots  to  separate  the  sequential  values  of  condi¬ 
tioning  envelope.  The  normalization  was  performed  relative  to  the  plot 
labeled  1 30  <  ^170  mv.  It  would  appear  from  these  data  that  the 

joint  channel  model  pdf  does  provide  a  useful  description  of  the  actual 
behavior,  in  particular,  the  model  predicts  the  sharp  reduction  in  slope 
of  the  density  tail,  compared  to  the  unconditioned  rf  noise  density  of 
the  previous  figures  in  this  chapter. 

We  have  noted  (section  3.  5.  1)  that  the  ratio  of  mean  to  rras  has 
often  been  used  to  characterize  the  non-gaussian  nature  of  atmospheric 
noise  envelope  waveforms.  This  ratio  for  a  guassian  distribution  is 
0.  797  cr,  while  smaller  values  would  indicate  a  pdf  tail  that  is  larger 
than  gaussian.  Expressions  (3.  26)  and  (3.  27)  give  the  variance  and 
mean  value  of  the  magnitude  of  the  conditional  Hall  noise  component  of 
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our  model,  including  the  effect  of  the  instrumentation  system's  dynamic 
range.  The  limiting  value  of  this  ratio,  for  large  values  of  ymax» 
is  seen  to  be 


Expression  (3.  27) 

Lim  . . . . . . . . . 

y  —oo  n/ Expression  (3.  26) 
Jmax  r 

v  «y 
o  ■'max 


For  the  entire  range  of  our  noise  model  (i.  e. ,  background  component 
to  Hall  component),  we  would  expect  that  the  ratio  of  the  mean  to  rms 
of  the  conditional  I  rf  I  noise  distribution  is  gaussian  for  small  values  of 
the  pilot  channel  envelope  where  only  the  gaussian  background  compo¬ 
nent  of  the  noise  is  present  (x(t)  =  0),  decreasing  to  the  lower  Hall  compo¬ 
nent  limit  given  above  and  then  rising  towards  unity  as  the  upper  limit 
ymax  is  approached  by  the  pilot  channel  envelope. 

The  Hall  component  ratio  was  evaluated  numerically  and  found  to  be 
relatively  insensitive  to  the  exact  values  of  and  tr^.  The  result  is 
shown  in  Figure  3-7  for  the  case  m  =  1,  along  with  several  plots  of  the 
actual  ratio  computed  from  the  joint  channel  data  for  m  =  1  noise  condi¬ 
tions.  We  see  that  the  data  ratio  follows  the  general  behavior  predicted 
by  the  model,  in  particular,  the  increase  towards  gaussianity  at  small 
values  of  the  pilot  channel  envelope  supports  the  concept  of  the  back¬ 
ground  gaussian  component  of  the  atmospheric  noise.  The  actual  point 
of  this  transition  is  dependent  upon  relative  channels  gains  and  hence 
has  no  significance.  The  subscript  at  each  rf  channel  frequency  is  the 
observation  bandwidth  in  kHz  and  we  see  that  there  is  a  definite  trend 
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to  lower  values  of  the  mean-rms  ratio  at  the  10  and  20  kHz  bandwidths. 
This  indicates  that  the  specific  form  of  the  joint  channel  pdf  is  less  effec¬ 
tive  at  these  bandwidths,  even  though  the  particular  data  match  used  in 
Figures  3-1  through  3-4  (conditional  <r  i  statistic)  agrees  quite  well. 

y  I  v0 


Figure  3-7.  Comparison  of  Mean-RMS  Ratio  for  Conditional 
RF  Noise  Magnitude  Density 


3.  7 _ N oise  Model  Time  Structure 

We  have  specified,  at  this  point,  the  time  structure  of  two  compo¬ 
nents  of  our  noise  model.  The  gaussian  processes  n^(t)  and  n2(t),  in 
the  model 


y(t)  =  n^t)  +  x(t)  a(t)  n2(t). 
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are  assumed  statistically  independent,  with  each  having  jointly  gaussian 
statistics  with  correlation  determined  completely  by  the  bandlimiting 
filter's  impulse  response.  Since  we  have  specified  x(t)  as  a  two-state 
process,  it  would  be  very  difficult  to  specify  a  time  structure  of  x(t)  a(t) 
that  would  generate  sample  records  exactly  like  the  type  seen  in  Fig¬ 
ure  2-4.  However,  the  fact  that  we  seek  principally  a  model  that  is  use¬ 
ful  in  radio  receiver  design  relieves  a  considerable  amount  of  the  problem 
of  exact  specification  of  short-time  structure.  The  reason  for  this  is 
that  all  VLF-LF  receivers  must  employ  sharp  bandlimiting  filters  to 
remove  adjacent  channel  interference.  Signals  designed  for  these  radio 
systems  must  provide  long  observation  times  relative  to  these  band- 
widths1^’  and  thus  a  noise  model  that  describes  the  average  be¬ 

havior  of  the  noise  waveform  over  these  time  intervals  is  adequate.  Our 
goal  in  this  section  is  to  postulate  a  time  structure  for  x(t)  a(t)  that 
approximates  the  observed  behavior  reported  in  Chapter  2  in  the  r  oder- 
ate  and  long-time  sense  and  does  not  compromise  the  first  order  model 
characteristics  given  in  3.6. 

From  the  time  structure  observations  reported  in  Chapter  2, 
we  recall  that  the  bandlimited  atmospheric  noise  waveform  has  a  power 
spectral  density  determined  by  the  bandlimiting  filter.  Measurements 
of  the  fluctuating  power  level  of  the  noise  (the  estimates  A(t.)  closely 
related  to  the  x(t)  a(t)  model  component)  showed  that  it  was  essentially 
uncorrelated  beyond  the  filter  correlation  time,  although  secondary  cor¬ 
relations  of  3-4  msec  and  300-500  msec  were  observed  and  related  to 
mechanisms  of  the  noise-generating  process.  In  terms  of  the  ultimate 
problem  of  receiver  design,  these  observations  indicate  that  the  first 
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order  pdf  of  the  noise  is  the  most  important  characteristic  to  be  mod¬ 
eled.  Where  the  fluctuating  noise  power  level  is  important  due  to  assumed 
signal  characteristics,  one  can  use  the  joint  channel  model  to  estimate 
the  instantaneous  value  to  x(t)  a(t)  more  accurately  than  an  attempt  to 
exploit  the  small  time  structure  component  in  the  A(tp  series.  Thus  our 
efforts  in  this  section  are  directed  at  completeness  of  the  model  and  the 
development  of  a  canonic  model,  including  time  structure,  suitable  for 
Monte  Carlo  computer  simulation,  rather  than  an  ultimate  use  of  the 
x(t)  a(t)  time  structure  in  the  receiver  design  problem. 

3_!_LlL  General  Time  Structure 

The  proposed  time  structure  of  the  x(t)  a(t)  component  of  our  noise 
model  is  the  following: 

a)  The  two-state  process  x(t)  is  Markov  with  transition  rates 

and  The  X.^  rate  is  constant  and  chosen  to  provide  the 

short-time  correlation  of  A(h).  The  X^^t)  rate  is  stochastic 
and  controls  the  time  intensity  of  the  non-gaussian  excursions 
caused  by  x(t)  =  1. 

b)  The  stochastic  rate  \gj(t)  is  driven  by  a  statistically  independent, 
two-state  Markov  process  w{t)  with  transition  rates  Pqj  and  p^. 
The  Pjq  transition  rate  is  constant  and  chosen  to  provide  the  long¬ 
time  correlation  associated  with  multiple  discharges.  The  Pqj 
rate  is  also  constant  and  is  chosen  to  provide  an  approximation 

to  the  magnitude  of  the  long-time  correlation  observed  in  the 
A{t.)  record. 

c)  The  process  a(t)  is  a  nonlinear  function  of  x(t):  a(t)  assumes  a 
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fixed  value,  described  by  the  random  variable  a  with  pdf  given 
by  (3.  2),  for  each  discrete  transition  of  x(t)  from  0  to  1.  Each 
succeeding  value  of  a(t)  is  statistically  independent  of  those  pre¬ 
ceding. 

We  see  that  the  effect  of  our  model  is  to  describe  each  burst  event, 
at  the  output  of  the  bandlimiting  observation  filter,  as  a  noise  burst  of 
constant  power  level  which  is  a  random  variable  for  different  burst  events. 
The  occurrence  intensity  of  these  discrete  noise  bursts  is  further  modu¬ 
lated  by  a  two-state  intensity  process,  approximating  the  occurrence  of 
single  discharge  events  and  the  near  continuum  of  noise  found  in  multiple 

discharges.  This  type  of  model  is  similar  to  that  originally  proposed  by 
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Furutsu  and  Ishida  and  recently  used  by  Coon  in  the  development 
of  an  analog  noise  simulator.  Their  approach  used  a  train  of  Poisson 
distributed  impulses  to  excite  the  bandlimiting  filter,  with  additional 
burst  packets  of  impulses  described  by  a  higher  intensity  Poisson  pro¬ 
cess  to  represent  multiple  discharges.  Our  noise  burst  model  has  a 
plausible  basis  in  the  leader  structure  of  single  discharge  events  which, 
at  frequencies  of  LF  and  above,  tend  to  appear  as  a  short  burst  of  noise. 
Where  the  burst  energy  is  more  time  concentrated,  we  are  simply 
approximating  the  effect  as  a  noise  burst  since  we  are  interested  in  the 
effects  on  a  scale  of  milliseconds  or  longer.  For  the  case  of  multiple 
discharges,  we  observed  in  Chapter  2  that  these  events  are  character¬ 
ized  by  an  additive  combination  of  a  continuum  noise  burst  and  shorter, 
high  intensity  bursts.  Our  model  tends  to  provide  this  behavior,  with 
the  continuum  noise  approximated  by  the  high  intensity  event  rate  and  the. 
random  variable  "a",  describing  each  discrete  power  level,  causing 
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large  excursions  to  model  high  intensity  bursts.  The  model  implicitly 
ignores  the  correlation  of  noise  burst  amplitude  with  time  intensity  as 
observed  for  multiple  discharges.  This  preserves  the  first  order  model 
as  presently  developed. 

3.  7.  2  Functional  Form  of  the  Model  Time  Structure 

The  assumed  form  of  the  time  structure  for  x(t)  a(t)  now  allows  us 
to  determine  the  functional  relationship  of  the  A(t^)  estimates  to  the 
model  and  hence  estimate  the  model  parameters  necessary  to  provide 
the  model  with  the  experimentally  observed  characteristics.  The  A(t.) 
estimates  were  generated  by  an  integrate -sample -dump  technique  which 
is  mathematically  equivalent  to  the  system  shown  in  Figure  3-8.  We  see 


Figure  3-8.  Equivalent  A(t^)  Generator 

from  the  well-known  relations  among  autocorrelation  functions  that  we 
can  determine  the  form  of  the  A(t.)  time  series  autocorrelation  by  cas¬ 
cading  the  effect  of  each  autocorrelation  function.  To  determine  the 
form  of  the  autocorrelation  of  v(t)  (where  we  neglect  double  frequency 


terms  from  the  linear  rectification  of  the  bandpass  y(t)  waveform)  we 
utilize  the  discrete  nature  of  the  Markov  x(t)  generator  to  write 


R  (t)  =  P(x,=l)  P.  (N  =0)  v|p  (t)  +  P  (x  =0)  P.  (N=0)  VpP  (t) 
v  z  *10  T  n  v  xt  *01  w  V 


+  Px  .  (N-even)  v£  +  P.  .  (N-odd)  vr 
*01,AT0  U  K01*K10  ° 


VH 


+  P(x,  =1)  P.  .  (N-even)  v„  +  P.  ,  (N-odd)  vn  vw  . 

t  N)1,X10  H  \)1'X10  G  H 

(3.30) 

In  the  first  term  of  (3.30)  P(xt=l)  is  the  probability  that  x(t)  =  1,  which 

is  stochastic  and  slowly  varying  as  determined  by  the  second  Markov 

generator  (3.7.  1 -(b)  above).  P.  (N  =0)  is  the  t  dependent  probability 

10  T 

that  zero  Poisson  distributed  transitions,  at  rate  occur  in  time  t, 

_ 

and  vH  Py(T)  is  the  autocorrelation  of  the  Hall  component  envelope.  The 

function  pv(t)  is  the  normalized  autocorrelation  of  a  linear  envelope 

34 

detector  driven  by  a  narrow-band  gaussian  process.  The  remaining 
terms  represent  the  other  possible  combinations  of  Poisson  events  at 
time  t  +  t  given  x(t)  =  1,0  at  time  t.  These  terms  may  be  regrouped 
in  order  of  importance  to  yield 


R  (t)  =  p  e 
v  *x 


'xioT 


vHpv<T)  +  fl(7H)  +  f2(vGvH> 


+  f3(4)  +  f4(vGe  10  ) 


(3.31) 


The  first  term  then  is  dependent  upon  the  variance  of  the  Hall  envelope 
component,  the  autocorrelation  of  the  noise  envelope  and  the  exponential 
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decay  of  the  P.  (N-0)  term.  The  second  term  represents  uncorrelated 
10 

averages  of  Vpj  caused  by  x(t)  transiting  an  even  number  of  times 
during  time  t,  while  the  last  terms  represent  much  smaller  effects. 

Neglecting  for  the  moment  the  time  dependence  of  px(t),  which  will 
be  used  to  model  the  300-500  msec  multiple  discharge  effects  and  hence 
will  not  affect  the  short  term  behavior,  and  assuming  the  Pv(t)  function 
decays  rapidly  to  a  Pv(°°)  value  (valid  for  1  kHz  and  larger  rf  bandwidths 
with  t  >  1  msec),  we  see  that  the  sampled  autocorrelation  values  are 
given  by 


RA(t)(nT) 


nT+T 


'nT-T 


px(t)  e  10  v*p  («o> 


Rv(t) 


1  -  |nT-T 


Rh(nT-T) 


dT 


+  constant  terms,  n  ^  1 


RA(t.)(nT)  =  px(t>  VHpv<”) 


-(nT+T)\ 


10 


(>r +  nT) 

\\i0  / 


-(nT-T)\;o  /  i  \  2  ’nTK10,  TX  .  n 

+  e  \K~n  +  nT)-\Ze  (nTX10+l) 


TO 


10 


+  constant  terms. 


(3.32) 


Thus  the  autocorrelation  of  the  model  A(t.)  series  is,  for  short-time, 
deoendent  upon  first  order  amplitude  statistics  previously  specified 
^ Px(t)  and  v^  and  an  exponential  decay,  with  parameter  Xi0  of  the 
Markov  x(t)  generator.  Physically,  the  parameter  \1Q  is  related  to  the 
average  length  of  single  discharge  events,  which  is  of  the  order  of 
1  msec.  The  exact  choice  of  this  parameter  will  be  based  on  simulation 
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results  given  in  section  3.  7.  3  below. 

For  correlation  times  nT  »  the  dominant  short-time  term  of 

So 

(3.  32)  will  decay  to  zero  and  the  A(t.)  autocorrelation  will  be  given  by 


_2 


RA(t  )(nT)  =  px(t)  px(t+nT)  v^  +  second  order  terms  (v^*  v^,  Vq) , 


(3.33) 


where  we  treat  R^t)  as  a  unit  delta  function  with  respect  to  the  slow 
variation  of  p(t)  p(t+nT)  which  will  be  of  the  order  of  hundreds  of  milli¬ 
seconds.  The  first  term  of  (3.  33)  is  the  first  term  of  the  second  square 

brackets  of  (3.  30)  where  we  have  used  the  nT  »  assumption  to 

X10 

neglect  transients  in  the  x(t)  generator.  The  joint  probability,  P(x(t)  = 
l,x(t+nT)  =  I),  will  only  be  influenced  by  the  stochastic  intensity  ^(t). 

To  evaluate  p  (t)  p  (t+nT)  we  make  the  following  definitions  and 
observations  concerning  the  Markov  generated  w(t)  which  controls  kQ1(t): 
a)  w(t)  is  the  state  of  the  second  Markov  generator  controlling 
the  x(t)  intensity.  w(t)  has  rate  parameters  |Xq1  and 


b)  w(t)  •!  —  **,  -  Px  =  J 


01 


\ii  +  No 


W(t)  =  0  ~~  -  p l  = 


01 


01  Fx  vs  .  . 

\>i +  No 


Hn 1 

c)  P  -  P(wf=l)  =  — - 7— - 

w  t  H01  +  F-10 


q  =  1  -  p  , 
w  *w 


k  =  H10  +  N)1  =  ^10^qw’ 


d)  The  w(t)  state  equations  are  given  by 
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p(w(t)=0) 

-k-r  , 

p  e  +  q 

rw  Mw 

qwd-e-kT) 

p(w(0)=0) 

p(w(r)=l) 

/t  -ki\ 

_  pw  "e  > 

-k-r  , 

qw  e  +  p 

W  w 

p(w(0)=l) 

1 


The  p  (t)  p  (t+nT)  term  with  nT  »  r —  ,  can  now  be  evaluated  using 
x  x  \10 

the  discrete  states  of  w(t)  as 


p(t)  p(t+nT)  =  P(xt.xt+nT=  1 1'Y'VnT11  0  P(,t'”t+nTs  0 


+  P(Y  xt+nT  =  1 1  wt  =  0,  w(+nT  =  1 )  P(w,  =  0,  wt+nT  =  1 ) 

+  p<xt,xt+nT'l|wt=>.wl+nT  =  0)- 
p(wt  =  1 .  wt+nT  •  0)  +  P(xt.  xt+nT  =  1 1  wt.  wt+nT  =  0)  P(w,.  wt+nT  =  0) . 

Using  the  chain  rule  for  conditional  Markov  probabilities  and  the  state 
equations  given  above,  the  expression  becomes 

Px<0  PxU+nT)  =  (p'  )2  pw(qw  e'knT  +  Pw)+  p'p^(pw-q„,  e_knT) 

+  pxPxpw (%•'>»  e"knT)+  (px)2  %(Pw  e‘knT  +  O' 

This  can  be  rearranged  to  yield 

px(t)  px(t+nT)  =  e‘knT  ((p')2  pw  +  (pk)2  qw  -  P^  )  +  P*.  «•  34) 


p  -  p  (t)  =  p  p^  +  q  pS. 
rx  *x  rwrx  wrx 


where 
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We  recognize  the  left  term  of  (3.34),  which  is  the  dynamic  portion,  to 

f  s 

be  the  variance  of  a  two  state  process  with  state  values  p^  and  px  and 
probabilities  pw  and  q  ,  respectively,  while  the  right  term  of  (3.  34)  is 
the  average  value  squared  of  the  process.  Combining  (3.  34)  with  (3.  33) 
we  ha'^e  the  long-time  behavior  of  the  A(t.)  mode]  autocorrelation  given 

by 


Ra(U(nT>  ■ 


-knT 


(\(t>)  +  ?x_ 


_2  2  _  _ 

v^  +  second  order  terms  ^g,vGvH^‘ 


(3.35) 


The  model  time  structure  provides  five  independent  parameters.  The 
first  two,  Xjq  and  k  =  K11q/qw  are  chosen  to  provide  the  approximate  ex¬ 
ponential  decay  of  the  data  for  short-time  and  long-time,  respectively. 
The  difference  in  the  form  of  these  two  decay  parameters  is  caused  by 
the  respective  sources  of  the  correlation.  For  short-time,  the  correla¬ 
tion  depends  upon  the  length  of  each  noise  burst  which  is  the  probability 
of  zero  Poisson  distributed  transitions  in  the  time  interval  nT,  whereas 
the  latter  depends  only  upon  x(t)  being  in  state  1  at  the  beginning  and  end 
of  the  time  interval,  independent  of  the  number  of  intervening  transitions 
or  "path"  of  x(t). 

The  remaining  three  parameters  are  related  by  two  equations  and 
an  inequality.  The  average  probability  of  x(t)  must  equal  that  required 
by  the  first-order  model, 

’  »model  *  ?w?x  +  Vx  =  V  <3-  36) 

The  range  of  the  x(t)  probability  can  be  determined  from  the  data  as  a 
normalized  ratio 


px(t)2/p-(t)2 


[(px)  pw+(px)  % 


RA,Data<10  msec)/RA,  Data'”’  *  RDATA-  <3'37> 


and  finally,  we  have  the  normal  probability  constraint 


0<px'PxSI- 


(3.38) 


In  selecting  a  parameter  set  to  satisfy  these  relations,  two  cases 

f 

arise.  The  first  occurs  when  px  <  1  and  a  non-unique  solution  exists. 

f  s 

We  can  arbitrarily  select  p  -  q  ,  which  leads  to  symmetric  p  and  p 

W  W  X  X 

values.  It  can  be  shown  that  this  choice  minimizes  the  difference, 
f  s 

px  -  Px,  for  a  given  RDATA  ratio.  The  parameters  for  this  case  are 
pw  =  qw  =  °*  5‘  px  =  px  +  Px^ RDATA  -  ; ,  Px  s  Px  ~  Px^RDATA--T. 


X  ‘  X  ‘  X 


(3.39) 


It  is  also  possible  to  arbitrarily  assign  px  at  a  large  value  to  repre¬ 
sent  the  near  continuum  noise  of  multiple  discharges,  or  find  that  (3.  39) 

f 

results  in  p  >  l  for  given  p  and  RDATA  values,  in  which  cases  the 

X  X 

g 

remaining  two  parameters,  p,  ,  p  ,  are  given  by  the  solutions  of 

Pw[(px)2  -  Px]  +  Pw[px  -  Sv£  +  RDATA  Px]  +  [px  -  Px]  ■  °- 


p  ■  p  p 

*X  { XrW 
1  -  p 

”  u> 


(3.40) 


The  equation  set  (3.40)  reduces  to  two  linear  equations  for  the  second 


f 

case  where  px  =  1 .  This  is  representative  of  tropical  and  frontal  con¬ 
ditions. 

We  noted  in  our  discussion  of  the  observed  A(t)  sample  record 
(section  2.  5.  2)  that  the  behavior  suggests  a  slowly  varying  continuous 
modulation  of  a(t)  combined  with  the  uncorrelated  fluctuations.  The  two- 
state  intensity  modulation  proposed  here  is  a  coarse  approximation  to 
this  and  one  might  wonder  why  not  modulate  the  intensity  of  x(t)  by  a 
continuous  lowpass  random  process.  We  explor  ed  this  idea  through 
computer  simulations  of  A(t^)  and  found  that  the  technique  works  well 
for  moderate  values  of  px#  For  the  larger  values  of  px,  .6  <  px  <  .9, 
as  required  for  tropical  and  frontal  conditions,  the  nonlinear  relation 
between  modulation  of  X.jg(t)  and  px(t)  introduced  by  the  relation  px(t)  = 
\)1^/^01^  +  X10^  maltes  ^  difficult  to  control  either  px  or  the  effective 
correlation  structure.  The  former  effect  then  compromises  the  first- 
order  noise  model  which  is  very  undesirable.  For  these  reasons  we 
have  concluded  that  a  two-state  or  Poisson -Poisson  model  is  the  better 
form,  short  of  incorporating  amplitude  correlation  with  the  time  noise 
burst  intensity.  Our  model  for  x(t)  a(t),  then,  yields  a  functional  behavior 
similar  to  that  found  in  the  experimental  observations  and  can  be  relat  ‘d 
to  the  physical  mechanisms  of  the  atmospheric  noise  sources.  Further, 
the  time  structure  parameters  can  be  chosen  in  a  manner  consistent 
with  the  first  order  model  parameters  and  these  choices  do  not  depend 
in  an  absolute  sense  upon  moment  statistics  of  the  first  order  model. 

In  the  next  section  we  shall  demonstrate,  via  a  computer  simulation, 
that  the  combined  first-order  and  time  structure  parameters  lead  to 
simulated  A(t.)  sample  records  and  autocorrelation  that  agree  well  with 


observations. 


3.7.  3  Time  Structure  Simulation  Results 

To  verify  that  the  model  time  structure  produces  the  type  of  behavior 
predicted  in  the  previous  section,  and  that  the  results  are  compatible 
with  the  first-order  model  parameters,  a  Monte  Carlo  simulation  was 
run  for  several  cases.  The  details  of  the  complete  atmospheric  noise 
simulation  generator  are  given  in  section  3.8  and  Figure  3-11.  For  our 
purposes  here,  the  x(t)  and  v(t)  two-state  generators  were  simulated 
with  uniform  random  generators  mapped  to  an  exponential  pdf  to  produce 
the  time  of  the  next  switching  event.  Each  A(t)  1  -msec  sample  was 
composed  of  the  sum  of  five  Rayleigh-distributed,  statistically  indepen¬ 
dent,  samples  plus  additional  Rayleigh -distributed  samples 

R(ax<rvH)  for  the  time  period  that  x(t)  =  1.  The  second  Rayleigh  param¬ 
eter,  a  ,  was  generated  for  each  x(t)  0-1  transition  from  a  x(m,  1)  dis- 
tribution,  generated  according  to  (3.  1).  All  of  the  parameters  in  these 
distributions  were  taken  from  the  first  first-order  model  parameters 
summarized  in  Table  3-1  in  section  3.8. 

A  simulated  sample  record  of  the  A(t^)  estimates  is  shown  in  Fig¬ 
ure  3-9  with  tropical  first-order  parameters  and  time  structure  param¬ 
eters  described  below.  Comparison  of  this  record  with  Figure  2-7,  the 
observed  A(t.)  time  series,  shows  a  reasonable  imitation,  although  the 
distinct  undulating  behavior  of  the  observed  data  is  not  present  in  the 
simulation. 

An  autocorrelation  of  the  simulated  A(t.)  estimates  is  shown  in 
Figure  3-10.  The  two  decay  parameters  were  empirically  selected  to 
produce  the  results  shown  there.  Equations  (3.40)  were  used  to  deter¬ 
mine  the  remaining  three  parameters  as 


Figure  3-9.  Simulated  A(t.)  Sample  Record 


Figure  3-10.  Autocorrelation  of  Simulated  A(t.)  Sample  Record 


RDATA  estimated  by  (3.  37)  =  1. 22 
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=  .99.  -  pw  -  0.66,  -  p®  =  0.27. 

Since  the  short-time  decay  parameter,  Xjq,  is  independent  of  all  other 

parameters  and  our  observations  indicated  that  this  decay  was  found  for 

all  noise  conditions,  the  empirical  value  of  850  Hz  is  similarly  valid  for 

all  noise  conditions.  The  long-time  decay  parameter,  Pjq,  must  be 

computed  for  each  parameter  set  as  pt j q  =  2qw,  where  k  value  of  2  Hz 

is  assumed  to  be  representative  of  all  multiple  discharge  phenomena. 

The  model  time  structure  produces  the  results  shown  in  Figure  3-10  for 

all  noise  conditions  except  frontal  where  the  large  value  of  p  ,  combined 

with  the  model's  lack  of  amplitude  correlation  with  xit)  intensity,  pre- 

o 

vents  a  sufficiently  large  value  of  to  match  the  observed  autocorre- 

*x 

lation.  All  of  the  time  structure  parameters  are  summarized  in  sec¬ 
tion  3.  8  below. 

3.  8  A  Canonic  Atmospheric  Noise  Generator 

The  combined  first-order  and  time  structure  model  for  bandlimited 
atmospheric  noise,  including  the  joint  channel  model,  can  be  conveniently 
realized  for  digital  computer  simulation.  The  noise  generator  is  shown 
in  Figure  3-11,  where  the  portion  above  the  dashed  line  realizes  the  time 
structure  and  a(t)  of  the  Hall  component,  wheieas  the  lower  realizes 
the  actual  envelope  or  rf  noise  samples.  The  upper  portion  is  iterated 
each  time  the  state  of  that  portion  changes,  as  determined  by  the  random 
event-time  samples,  TNXT.  The  nonlinear  mappings  from  Random 
Number  Generators  (RNG's)  with  uniform  distribution  fU(0,  1))  are 
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conventional  inversion  formula  which  can  be  used  when  a  desired  pdf  has 

an  analytic  cumulative  distribution  from  which  the  inverse  mapping  can 
35 

be  obtained.  Other  generator  forms,  such  as  the  in-phase  and  quadra¬ 
ture  phase  components  of  the  bandlimited  atmospheric  noise  waveform 
can  be  similarly  constructed.  If  one  is  interested  only  in  statistically 
independent  atmospheric  noise  samples,  without  consideration  of  time 
structure,  the  Hall  component  of  either  envelope  or  rf  noise  or  the  joint 
samples  of  both  can  be  realized  with  U(0, 1)  RNG's  and  the  appropriate 
mapping.  These  mappings  can  be  determined  for  the  integer  values  of 
m  that  we  have  considered  in  our  first-order  model.  Figure  3-11  does 
not  show  the  truncation  of  sample  values  at  the  maximum  upper  limit  of 

Table  3-1.  Canonic  Noise  Generator  Parameters 
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the  various  pdf's.  This  can  be  accomplished  either  directly  as  the  Hall 
component  output  by  iterating  for  a  new  value  if  the  maximum  is  exceeded 
or  by  limiting  the  a.  generator  range  as  suggested  by  Hall.  While  these 
produce  slightly  different  pdf's  at  the  maximum  range  the  differences 
wouxd  probably  not  be  important  in  most  simulation  applications. 

The  parameters  required  by  the  model  and  canonic  generator  are 
given  in  normalized  form  in  Table  3-1.  This  table  is  a  representative 
summary  of  all  of  the  sample  records  that  we  analyzed  from  the  data 
collection  program  in  New  Jersey. 

3.  9  Summary 

3.  9. 1  Model  Results 

In  Chapter  1  we  briefly  described  the  nature  of  electromagnetic 
radiation  from  lightning  discharges,  which  is  the  principal  source  of 
low  frequency  atmospheric  radio  noise.  The  complex  time  structure 
of  these  discharge  radiations,  especially  at  frequencies  above  40-50  kHz, 
provide  motivation  for  a  multiplicative  model  of  the  noise.  y{t)  =  A(t)  n(t), 
as  originally  suggested  by  hall.  This  model  describes  the  noise  wave¬ 
forms  observed  at  the  output  of  a  bandlimiting  filter,  rather  than  employ¬ 
ing  the  more  difficult  approach  of  statistically  modeling  the  discharge 
process  itself  and  then  determining  the  resulting  filter  output  statistics 
from  this  source  model.  We  extended  Hall's  model  by  suggesting  that 
A(t)  could  be  considered  to  be  independent  of  frequency  and  thus  provide 
a  statistical  link  between  noise  waveforms  observed  in  different  frequency 
channels. 

In  Cahpter  2  we  reported  various  experimental  observations  made 
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at  VLF  and  LF  to  explore  what  statistical  and  dynamic  behavior  of  A(t) 
was  required  to  approximate  actual  noise  waveforms.  Chapter  3  has 
developed  two  aspects  of  a  mathematical  model  based  on  the  multiplica¬ 
tive  concept  and  compared  them  to  the  observed  data  to  determine  appro¬ 
priate  parameters.  Our  proposed  model  for  bandlimited  atmospheric 
noise  waveforms  is  given  by  y(t)  =  n^t)  +  a(x(t))  n2(t),  where  nj(t)  is  a 
statistically  independent  background  gaussian  process,  x(t)  is  a  two-state 
Markov  process  with  a  time -varying  transition  rate  parameter  and 
a(x(t))  n2(t)  is  the  multiplicative  process  given  by  Hall  where  the  x(t) 
process  acts  to  "turn-on"  or  "turn-off"  this  Hall  component  in  a  partic¬ 
ular  manner.  In  addition,  a(x(t»  is  independent  of  frequency  over  signif¬ 
icantly  large  frequency  increments  and  is  thus  identically  the  same  for 
noise  waveforms  observed  over  this  increment.  This  noise  model  repre¬ 
sents  a  synthesis  of  three  previously  suggested  models,  including  Hall's, 
with  the  addition  of  our  concept  of  the  frequency  independence  of  a(x(t)). 

Our  principal  comparison  of  this  noise  model  to  experimental  obser¬ 
vations  was  based  on  records  of  effectively  simultaneous  samples  of  the 
noise  envelope  in  one  channel  and  the  rf  noise  in  a  second,  both  in  the 
LF  band.  We  found  that  the  model  and  associated  six  parameters  would 
satisfactorily  describe  each  channel's  unconditional  probability  density 
and  a  conditional  statistic's  variation  linking  the  channels.  These 
six  parameters  compare  to  eleven  parameters  required  to  simply  fit 
mathematical  expressions  to  all  three  curves.  We  further  demonstrated 
that  the  joint  channel  model  (the  fact  that  a(x(t))  is  common)- predicted 
the  basic  form  and  relation  of  the  conditional  rf  noise  pdf  measured 
with  a  different  technique,  and  predicted  the  basic  relationship  of  a 


different  joint  channel  statistic  than  used  in  comparing  the  model  with 
the  data. 

The  basic  time  structure  of  our  model,  introduced  to  model  the 
observed  dynamics  of  A(t),  is  provided  by  the  stochastic  nature  of 
the  x(t)  transition  which  acts  to  "turn-on1’  the  a(x(t»  n2(t)  noise  burst. 
The  average  length  of  each  burst  is  appr  oximately  equal  to  that  reported 
for  individual  lightning  discharge  events.  This  behavior  of  a(x(t))  pro¬ 
vides  a  basically  uncorrelated  structure  for  A(t),  with  a  small  second- 
order  correlation  extending  over  several  milliseconds.  The  long-time 
correlation  observed  of  A(t)  is  approximated  in  the  model  by  modulating 
the  intensity  (Poisson  rate  parameter)  of  x(t)  with  a  second  two-state 
process  which  corresponds  roughly  with  long  multiple  discharge  phenom¬ 
ena.  The  model  also  provides  the  correct  autocorrelation  of  y(t)  itself 
via  the  time  structure  of  n^t)  and  n2(t),  although  this  fact  is  relatively 
unimportant  in  typical  low  frequency  signal  processing  design  problems. 

Finally,  we  interpreted  the  noise  model  relations  in  terms  of  a  block 
diagram  of  a  canonic  noise  generator,  suitable  for  Monte  Carlo  com¬ 
puter  simulations.  This  generator,  together  with  both  first-order  and 
time  structure  parameters,  was  used  to  demonstrate  that  a  simulated 
A(t.)  time  series  yielded  the  same  autocorrelation  as  that  of  the 
data. 

3,  9.  2  Additional  Questions 

The  experimental  program  described  here  was,  of  necessity,  an 
attempt  to  collect  specific  data  based  only  on  the  general  concept  of  a 
multiplicative  noise  model.  The  final  form  of  the  model  in  turn  suggests 
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a  number  of  more  specific  experimental  questions.  In  general  these 
relate  to  the  variation  of  the  model  parameters  with  band  center  fre¬ 
quency,  joint  channel  frequency  separation  and  observation  filter  band¬ 
width.  We  suspect  that  the  form  of  the  model  is  applicable  to  much 
higher  frequencies  than  those  tested  here,  principally  because  the  leader 
structure  of  the  lightning  discharges,  which  appear  as  a  noise  burst, 
becomes  dominant  with  increasing  frequency.  Appropriate  questions 
regarding  joint  channel  separation  might  involve  more  precise  measure¬ 
ment  of  conditional  densities  at  varying  separations,  using  matched 
filter  envelope  responses,  and  perhaps  experimentation  with  other  chan¬ 
nel  waveforms  or  the  use  of  more  than  one  additional  channel. 

Questions  relating  to  the  effect  of  the  observation  filter  bandwidth 
are  perhaps  the  most  important  since  this  enters  intimately  into  any  radio 
system  design.  We  have  used  bandwidths  of  1,  10,  and  20  kHz  and  have 
not  found  significant  differences  in  the  first-order  model  parameters  at 
these  bandwidths,  although  the  results  of  section  3.6.  3  indicated  that  the 
specific  form  of  the  conditional  rf  noise  density  seemed  less  appropriate 
at  the  larger  bandwidths.  A  check  of  the  absolute  values  of  and  crnH 
used  in  the  tropical  model  at  1  and  10  kHz  showed  that  they  differed  by 
approximately  the  ratio  of  the  square  root  of  the  noise  power  bandwidths 
of  the  two  filters  used.  This  is  the  result  one  would  expect  for  a  process 
that  appeared  as  bursts  of  noise,  relative  to  the  filter  impulse  response. 
Conversely,  this  would  imply  that  one  could  scale  the  model  parameters 
given  in  Table  3-1  for  differing  bandwidths,  using  this  bandwidth  factor, 
at  least  above  1  kHz.  However,  as  the  bandwidth  is  reduced  significantly 
below  this  value,  the  time  width  of  individual  discharges  will  approach 
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the  time  constant  of  the  filter  and  clearly  this  will  affect  the  model 
parameters  in  a  nonlinear  manner  since  the  lower  limit  of  such  band¬ 
width  reduction  must  be  a  gaussian  noise  process.  In  terms  of  the  quali¬ 
tative  changes  in  our  model  parameters,  we  would  expect  that  the  rela¬ 
tive  level  of  the  gaussian  background  process  would  increase  and  the 
value  of  px  would  decrease  with  decreasing  rf  bandwidth.  This  means 
that  a  noise  pdf  such  as  associated  with  tropical  conditions  would 
trend  towards  that  found  for  quiet-night  or  quiet  with  a  decreasing  rf 
bandwidth. 

This  type  of  change  in  noise  pdf  is  important  in  signal  processing 
design  for,  as  will  be  seen  in  the  next  chapter,  there  is  an  8  dB  differ¬ 
ence  in  signal-to-noise  ratio  improvement  that  can  be  obtained,  with  a 
simple  nonlinearity,  for  tropical  as  compared  to  quiet  noise  conditions. 
Intuitively,  the  wider  rf  bandwidth  preserves  the  effects  of  individual 
discharges  whereas  the  smaller  bandwidth  tends  to  blend  the  lower  level 
effects  of  these  together,  contributing  to  the  rise  in  the  apparent  back  • 
ground  level.  In  terms  of  signal  processing,  the  lower  level  blending 
"gaussianizes"  the  noise  in  those  time  regions  and  prevents  the  pro¬ 
cessing  structure  from  discriminating  against  the  effects  of  the  original 
discharges.  In  general,  then,  one  would  want  to  use  the  smallest  rf 
bandwidth  at  which  significant  changes  in  the  noise  pdf's  begin,  and  our 
observations  do  not  provide  this  information.  However,  since  there 
are  relatively  few  changes  in  pdf  character  from  I  to  10  kHz,  our  noise 
model  and  supporting  parameter  table  does  represent  the  "best"  (in  a 
signal  processing  sense)  noise  characteristics  one  could  obtain  from 
such  an  optimum  bandwidth  selection. 
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In  the  next  chapter  we  shall  apply  the  noise  model  to  the  design  and 
performance  analysis  of  a  phase -locked -loop  for  use  in  low  frequency 
navigation  systems.  The  predicted  performance,  based  on  the  model, 
will  be  verified  using  computer  simulations  and  sample  records  from 
various  noise  conditions. 


Chapter  4 

AN  OPTIMUM  DESIGN  FOR  NAVIGATION  RECEIVERS 
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One  of  the  most  significant  problems  in  the  design  of  low  frequency 
radio  receivers  is  the  development  of  a  performance  analysis  and  opt’ mi  - 
zation  which  includes  the  effects  of  non-gaussian  atmospheric  noise. 

This  development  requires  a  quantitative  description  of  noise  character¬ 
istics.  We  have  endeavored  in  the  preceding  three  chapters  to  develop 
such  a  quantitative  model,  with  associated  parameter  sets,  to  describe 
certain  aspects  of  bandlimited  atmospheric  noise.  We  appealed  at  sev¬ 
eral  points  in  that  work  to  our  ultimate  goal  of  developing  a  noise  model 
that  was  useful  in  the  receiver  design  problem.  This,  for  example, 
justified  our  focus  on  noise  characteristics  at  the  output  of  the  bandlimiting 
filter  and  provided  our  motivation  for  development  of  the  joint  channel 
noise  model.  In  this  chapter  we  shall  use  this  model  as  a  tool  to  analyze 
several  typical  navigation  receiver  processors  and  to  specify  a  near 
optimum  receiver. 

We  shall  assume  a  linear  signal  generator  model  and  use  a  sampled 
data  phase -locked -loop  (PLL)  as  the  basic  navigation  receiver  processor 
structure, with  a  linear  time -invariant  loop  providing  a  performance 
reference  for  our  analysis.  This  loop  can  be  optimized  for  the  atmo¬ 
spheric  noise  power  level  and  for  a  statistically  independent,  gaussian 
distributed  noise  sequence,  such  a  loop  would  be  optimum  in  the  mean- 
square  error  sense  without  restriction  as  to  structure.  We  have  seen 
that  a  principal  characteristic  of  atmospheric  noise  is  that  its  first- 
order  probability  density  is  not  gaussian.  To  account  for  this  we  shall 
introduce  a  zero  memory  nonlinearity  in  the  loop  and  show  that  we  can 
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determine  the  approximate  performance  of  such  a  loop  by  using  only  the 
first-order  probability  density  of  the  noise  model.  We  shall  use  the  noise 
model  pdf  to  compare  the  performance  of  several  nonlinearities  such  as 
a  hard  limiter,  clipper,  and  hole  puncher  to  the  optimum  performance 
achievable  with  this  loop  structure. 

Our  experimental  observations  have  shown  that  not  only  is  the  first- 
order  probability  density  of  atmospheric  noise  non-gaussian,  but  that 
samples  of  the  noise,  separated  by  a  spacing  greater  than  the  bandlimiting 
filter  correlation  time,  are  uncorrelated  but  are  not  statistically  inde¬ 
pendent  (this  produces  the  correlation  of  the  A(t^)  sequence,  in  some 
cases,  of  0.  5  to  1  second).  To  account  for  this  in  the  navigation  pro¬ 
cessor  we  shall  use  the  joint  channel  noise  model  and  the  simultaneous 
pilot  channel  envelope  observations  to  approximate  the  atmospheric  noise 
as  a  gaussian  process  with  a  time-varying,  random,  power  level.  The 
filter  which  is  unconditionally  optimum  in  the  mean-square  error  sense 
will  then  be  given  by  a  time-varying  Kalman-Bucy  filter.  This  filter  will 
optimally  compensate  for  the  non-gaussian  probability  density  of  the 
atmospheric  noise  and  for  the  noise  time  structure  and  its  interaction 
with  signal  dynamics.  Finally,  all  of  these  filters  will  be  simulated  on 
a  computer  and,  for  two  types  of  simulated  stochastic  signal  processes, 
the  filter  performance  will  be  determined  by  using  recorded  noise  se¬ 
quences  for  all  types  of  noise  conditions,  from  nearly  gaussian  to  the 
most  severe  noise  associated  with  frontal  weather  systems.  These  tests 
will  show  that  the  simple  nonlinear  loop,  with  a  practical  nonlinear 
element,  provides  performance  improvements  of  6  to  16  dB  compared 
to  the  time-invariant  linear  loop,  and  that  for  loop  time  constants  less 
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than  several  seconds  the  near-optimum  time -varying  filter  provides 
further  improvements  of  3  to  5  dB  over  the  simple  nonlinear  loop. 

4.  1  Definition  of  the  Problem 

4. 1 . 1  Navigation  Problem 

Phase  coherent  navigation  systems  transmit  radio  signals  from 
two  or  more  known  geographic  locations,  where  all  of  the  system  signals 
are  at  least  phase  stable  with  respect  to  each  other  and  today  are  often 
synchronized  to  a  universally  accepted  time -frequency  scale  through  the 
use  of  atomic  resonance  oscillators.  A  receiver  on  a  mobile  platform 
receives  a  signal  from  each  station  which  is  given  by 

r(t)  =  y(t)  +  B(t)  sin  (u>ct  +  f*0a<  x(t),  V©  >  dT  +  Zero  Reference  Phase), 

where  y(t)  is  the  additive  atmospheric  noise,  *he  vector  inner  product 
yields  a  scaled  platform  velocity  component  perpendicular  to  lines  of 
constant  transmitted  phase,  and  the  "Zero  Reference  Phase"  represents 
a  reference  geographic  position  in  phase  coordinates.  The  two  types  of 
signal  formats  currently  in  operational  use  are  time  multiplexed  (LORAN) 
in  which  B(t)  is  a  bandlimited  pulse  waveform  which  each  station  trans¬ 
mits  at  a  prescribed  sequential  time  at  the  same  carrier  frequency,  or 
frequency  multiplex  (DECCA,  OMEGA)  in  which  B(t)  is  a  constant  for 
each  station,  which  transmits  at  a  different  frequency.  The  engineering 
of  a  complete  receiver  system  of  a  given  type  is  a  complex  problem 
involving  questions  such  as  determination  of  the  "Zero  Reference  Phase," 
acquisition,  interference,  noise  reduction  and  the  slructure  of  an 
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optimum  geographic  position  estimator  from  received  phase  estimates. 
We  shall  be  concerned  only  with  the  problem  of  estimating  the  received 
phase  from  one  station  when  the  received  carrier  is  corrupted  by  addi¬ 
tive  low  frequency  atmospheric  noise. 

4.  1. 2  Signal  Generator  Model 

To  model  the  dynamics  of  the  received  phase  component  perpendic¬ 
ular  to  lines  of  constant  transmitted  phase,  we  shall  assume  that  this 
may  be  prepresented  by  a  linear  mapping  from  a  vector  Markov  process 
as  follows: 

a)  Random  phase  process  generator  (s(t)  in  radians) 

s(t)  =  F  s(t)  +  G  u(t),  E[u(t)  u(t+T)]  =  Q6(t-T)  radian2,  E[u(t)]  =  0, 

b)  Received  phase, 

Sj(t)  =  [l  0  .  .0]  s(t)  =  H  s,  H  volts/ radian, 

c)  Atmospheric  noise  process,  y(t), 

E[y(t)  y(t+T)]  =  R(t-T)  volts2,  E[y(t)]  =  0, 

d)  Received  carrier  and  signal-to-noise  ratio, 

B(t)  s  l,  SNR  =  |r 

The  matrices  F  and  G  can  be  chosen  as  a  model  of  platform  dynamics 
plus  oscillator  dynamics.  We  shall  consider  only  sampled  data  systems 
in  which  the  sampling  interval,  AT, is  equal  to  or  larger  than. the  corre¬ 
lation  time  of  the  bandlimited  atmospheric  noise  process,  whence 
R{t-r)  =  R6(t--r).  As  we  noted  in  section  3.  7,  this  is  not  a  restrictive 


-118- 

assumption  since  all  low  frequency  receivers  must  employ  bandlimiting 
to  remove  adjacent  channel  interference.  Thus  the  performance  of  con¬ 
tinuous  systems  can  be  approximated  by  an  equivalent  sampled  system 
with  AT  set  equal  to  the  noise  correlation  time.  For  the  case  of  LORAN, 
a  time  multiplexed  system,  a  sampled  data  model  is  an  exact  descrip¬ 
tion  of  the  receiver  processing  that  must  be  employed.  Using  the  state 
transition  matrix  of  the  continuoi  s  time  signal  model,  4>(t,  tQ),  the  equiv- 
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alent  sampled  data  signal  model  is  given  by 

s(k+l)  =  $(AT,0)  sk+Ju(k)  (4.1) 

Elu(k)  u(k)T]  =  Q  J0AT  4(AT,  t)  G  GT  *(AT,  t)T  dr  (4.  2) 

r(k+l)  =  H  s(k+l)  +  y(k+l)  (4.3) 

y(k)  —  E[y(k)  y(k+n)]  =  R6(k-n)  =  2s?Jr  6^k"n^*  (4.4) 


Specification  of  the  continuous  time  signal  process  matrices  F  and 
G  (assumed  time -invariant)  then  determines  the  equivalent  linear  differ¬ 
ence  equations  (4. 1)  ior  the  sampled  data  model  and  the  covariance 
matrix  (4.  2)  of  the  driving  random  process.  We  note  in  equation  (4.4) 
that  we  have  parametrized  the  noise  power  in  terms  of  received  signal - 
to  -noise  power  ratio,  where  the  signal  is  the  rf  carrier,  not  the  phase 
message  process  s(t)  which  in  general  is  nonstationary.  This  is  the 
sampled  signal-to-noise  ratio  that  one  would  measure  at  the  output  of 
the  antenna  bandlimiting  filter. 


4.  1. 3  Signal  Processing  Structure 
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In  considering  signal  processing  structures,  we  shall  assume  that 
the  error  in  the  local  estimate  of  the  received  phase  message  is  small. 
Under  this  condition  it  is  well  known  that  a  phase-locked -loop  (FLL) 
provides  a  statistically  optimum  estimate  for  the  case  of  the  corrupting 
noise  being  white  and  gaussian.  Most  important  for  our  purposes  is 
the  fact  that  such  a  PLL  in  sampled  data  form  can  be  directly  interpreted, 
via  state  variables,  as  a  Kalman-Bucy  linear  estimator  which  brings 
considerable  insight  and  useful  results  to  bear  on  our  problem.  The 
majority  of  the  research  into  PLL  performance  has  been  directed  to¬ 
wards  analysis  when  the  linearized  assumption  does  not  hold  and  prob¬ 
lems  such  as  acquisition,  cycle  slipping  and  threshold  behavior  become 
the  dominant  considerations.  In  the  navigation  context,  however,  these 
problems  (except  acquisition)  are  largely  irrelevant  because  if  cycle 
slipping  has  any  significant  probability  of  occurring  the  phase  estimate 
becomes  useless  as  a  measure  of  position  since  one  does  not  know  which 
is  the  correct  rf  cyclefor  the  Zero  Reference  Phase  ,thus  losing  the 
measure  of  absolute  position.  Our  linearized  restriction  is,  in  fact, 
the  only  case  of  practical  interest  in  navigation  receiver  design.  Our 
problem,  for  low  frequency  systems,  is  that  the  assumption  of  gaussian 
noise  does  not  hold,  as  we  have  demonstrated  in  Chapters  1,  2,  and  3, 
and  we  require  design  and  analysis  techniques  to  specify  an  improved 
or  optimum  signal  processing  structure  for  these  atmospheric  noise 
conditions. 

For  a  sampled  data  system,  the  mixing  operation  of  a  PLL  in  which 
the  local  reference  is  multiplied  by  the  received  signal  to  generate  the 
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lowpass  phase  error  is  performed  by  sampling  the  received  signal  at 
tee  zero  crossing  time  of  the  local  reference  oscillator  whose  phase  is 
the  estimate.  If  our  local  phase  estimate  is  given  by  &(t),  the  sampling 
times  are  defined  by  solutions  of 

2(k+l)  7T  -  s,(k+l  |k) 

sin(w  t+s(t»  =  0  -  t  =  - - -  , 

c  «c 

at  which  times  the  received  signal  samples  are  given  by  (using  the  small 
error  assumption) 

r(k+l)  =  y(k+l)  +  sin  (2(k+l)  v  +  s^k+l)  -  s^k+l  | k)) 
a  y(k+l)  +  Sj(k+1)  -  Sj(k+1  |k) 

We  note  that  this  is  also  the  sampled  loop  error  and  we  introduce  the 
following  notational  definitions: 

a)  |(k+l  |k);  local  estimate  at  k+1  using  observations  to  k 

b)  s(k+l  jk)  =  ^(k+1)  -  Sj(k+1  |  k) ;  local  estimate  error  at  k+1 
using  observations  to  k 

c)  z(k+l)  =  y(k+l)  +  s(k+l|k),  actual  error  signal  at  k+1  using 
observations  to  k 

d)  s(k+l  |k+l),  local  error  at  k+1  using  k+1  observations. 

The  design  problem  now  is  to  determine  the  transformation  from  the 
sequence  of  loop  error  samples,  z(k+l)  to  the  estimates  j>(k+l  |k+l). 

We  know  that  if  the  noise  sequence  is  uncorrelated,  then  the  linear  time- 
invariant  minimum  mean-square  error  estimator  is  given  by  the  steady 
state  Kalman-Bucy  filter  in  v/hich  the  error  process  is  weighted  and 
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applied  to  a  model  of  the  signal  generator  as  a  correction.  The  signal 
generator  model  state  then  provides  the  updated  local  estimate 
s(k+l  |k+l),  which  in  turn  is  mapped  to  the  local  estimate  at  the  next 
sampling  time.  These  relations  are  summarized  in  Figure  4-1,  which 
includes  details  of  the  simulation  to  be  discussed  in  sections  4.  5  and 
4.6  below. 

We  see  in  Figure  4-1  that  the  loop  error  process  passes  through  a 
block  labeled  ZNL,  Zero  Memory  Non-Linearity,  and  a  possibly  time- 
varying  gain  matrix  K(k+1).  In  the  following  sections  we  shall  consider 
the  rationale  for  using  the  ZNL  at  this  location  and  develop  an  approxi¬ 
mate  method  of  determining  the  complete  estimator  performance.  The 
time-varying  gain  matrix  K  will  also  be  driven  by  a  real-time  computa¬ 
tion  of  the  matrix  variance  equation  based  on  using  the  pilot  channel 
envelope  sample  as  an  estimate  of  the  time -varying  power  level  of  the 
noise  sequence  y(k+l).  This  estimator  structure  provides  a  near  opti¬ 
mum  estimate  in  the  minimum  mean -square  error  sense  (MMSE)  with¬ 
out  restriction  as  to  estimator  class. 

4.2  Optimum  Linear  Time -Invariant  Estimator 
4.2.1  Kalman-Bucy  Equations 

The  experimental  observations  of  atmospheric  noise  waveforms 
reported  in  section  2.  5.  1  indicate  that  the  noise  process  is  uncorrelated 
beyond  the  correlation  time  of  the  bandlimiting  filter.  Samples  of  this 
process,  with  the  sampling  interval  meeting  this  requirement,  are 
similarly  uncorrelated.  Thus,  with  the  specified  signal  process  model 
and  this  observation,  the  MMSE  linear  filter  is  known  to  be  the 


Figure  4-1.  Block  Diagram  of  Phase -Locked  Loop  Simulation 
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Kalman-Bucy  filter.  For  a  gaussian  noise  process  and  our  assumed 
signal  model,  the  Kalman-Bucy  (K-B)  filter  is  the  MMSE  estimator  of 
all  possible  types  of  filters.  For  our  non-gaussian  atmospheric  noise 
process,  the  first-order  amplitude  statistics  are  non-gaussian  and  sam¬ 
ples  spaced  beyond  the  filter  correlation  time  are  uncorrelated  but  not 
statistically  independent  (our  model  attempts  to  describe  the  power  level 
correlation  of  such  samples).  Hence  the  K-B  filter  is,  in  our  case,  only 
an  optimum  filter  of  the  linear  class  although  it  is  easy  to  implement, 

t 

requiring  only  the  signal  model  and  unco  elated  noise  assumptions  plus 
a  measure  of  the  power  level  of  the  atmospheric  noise.  Sections  4.  3  and 
4.  4  to  follow  will  attempt  to  exploit  the  two  differences  noted  above  for 
atmospheric  noise  to  improve  upon  the  linear  filter,  although  we  will  find 
that  design  complexity  increases  accordingly. 

The  Kalman-Bucy  filter  equations  are  listed  here  to  establish  nota- 
tional  conventions. 

a)  Continuous  Time  Equations 

i)  Estimator 

4(t)  =  F  s(t)  +  K(t)[r(t)  -  H  s(t>]  (4.  5) 

ii)  Gain 

K(t)  =  S(t)  HT  R_1(t)  (4.6) 

iii)  Error  Covariance 

S(t)  =  F  2(t)  +  Z(t)  FT  -  2(t)  HT  R_1(t)  H  S(t)  +  G  Q  GT  (4.  7) 


iv)  Steady  State 
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K  =  S(oo)HTR_1  (4.8) 

where  S(oo)  is  a  positive  definite  solution  of 

0  "  F  S(oo)  +  s(oo)  FT  -  S(oo)  HT  R"1  H  S(oo)  +  GQG  (4.9) 

b)  Sampled  Data  Equations 

i)  Estimator 

(®  -  $(AT,  0))  (4.10) 

f(k+l  |  k+1)  =  S  s(k|  k)  +  K(k+1)[  z(k+l)]  (4. 1 1) 

ii)  Gain 

K(k+1)  =  £(k+l|k+l)  HT  R*(k+1)  (4.  12) 

iii)  Error  Covariance 

E(k+l|k)  =  *  S(kjk)  4T  +  GQ(AT)  GT  (4.13) 

S(k+1 1 k+1)  =  2(k+l  |k)  -S(k+1  |k)HT[R(k+l)  +  H  £(k+l  | k)  HT]_1  H  S(k+1  |k) 

(4. 14) 

iv)  Steady  State 

K  =  S(oo)HTR-1  (4.15) 

where  2(oo)  is  a  positive  definite  solution  of 

S(co)"1  =  HT  R-1  H  +  [$  E(oo)  $T  +  G  Q(AT)  GT]-1 


4.  2.  2  Sampled  Data/Continuous  Time  Relations 


-122- 


We  have  included  the  continuous  time  equations  above  because  it  is 
often  easier  to  solve  these  steady  state  equations  than  the  sampled  data 
equations.  For  reasonably  dense  sampling  with  respect  to  the  estimator 
filter's  impulse  response  there  should,  of  course,  be  a  correspondence 
between  these  two  solutions.  For  the  signal  model,  we  have  defined  the 
sampled  equations  in  terms  of  continuous  time  parameters  and  hence 
the  correspondence  is  immediately  established  via  F,  G  and  Q.  For 
the  white  observation  noise  samples,  we  have  a  time  series  with  a  known 
variance,  however,  the  variance  of  a  continuous  time  white  noise  pro¬ 
cess  is  infinite.  The  difficulty  is  that  the  definition  of  a  white  noise 
process  in  terms  of  an  impulse  correlation  function  is  only  useful  in 
terms  of  that  .unction's  action  on  an  integral  operator.  Thus  we  must 
seek  an  equivalent  white  noise  definition  that  results  in  the  specified 
noise  sequence  variance  when  acting  over  the  sampling  interval  AT. 

This  can  be  done  as  follows: 

y(t)  —  E[y(t)  y(t-Hr)]  S  R  6(t) 

i  rAT 

y(k)  =  at  J0  y(r)  dT 


ATE. 


•  • 


(4. 16) 
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This  relation  for  the  equivalent  area  of  the  continuous  time  white  noise 
correlation  function  will  be  shown  in  section  4.6. 1  to  yield  results  that 
are  accurate  to  better  than  1%  in  predicting  the  steady-state  behavior 
of  sampled  data  equations  from  the  continuous  time  steady- state  solu¬ 
tions  when  AT  is  less  than  1%  of  the  estimator  filter's  time  constant. 

4.  2.  3  Relation  to  Classical  Design  Procedures 

39 

Classical  design  procedures  for  linear  navigation  PLL's  have 
generally  been  based  on  a  specification  of  a  Type  I  or  Type  II  servo  loop 
(single  or  double  integrator,  respectively)  where  the  loop  parameters 
have  been  chosen  to  meet  a  transient  response  specification.  Noise  per¬ 
formance  has  teen  evaluated  via  the  noise  power  bandwidth  of  the  result¬ 
ing  linear  transfer  function,  assuming  a  white  noise  input.  Our  results, 
when  a  single  or  double  integrator  signal  generator  model  is  specified, 
can  be  related  to  this  point  of  view  by  relating  a  given  PLL  bandwidth 
specification  to  an  equivalent  noise  model.  This  relationship  will  be 
pointed  out  below  when  the  two  signal  models  are  specified.  The  prin¬ 
cipal  advantage  of  our  approach  is  that  we  bring  to  bear  on  the  problem 
many  useful  results  and  insights  from  modern  communication  theory 
and,  at  the  same  time,  prescribe  an  unambiguous  method  of  testing  the 
performance  of  different  receiver  structures.  This  becomes  important 
when  we  consider  nonlinear  or  time-varying  PLL's  where  there  is  no 
general  relationship  between  a  deterministic  signal  response,  such  as 
an  impulse  response,  and  the  performance  of  the  loop  when  stimulated 
by  signal  plus  noise.  Our  testing  procedure  will  consist  of  measuring 
the  error  of  various  PLL  structures  in  tracking  a  simulated  stochastic 
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signal  in  the  presence  of  the  recorded  atmospheric  noise  sample  se- 
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quences.  In  a  classical  context  this  parallels  the  technique  of  mea¬ 
suring  the  small  signal  loop  bandwidth  by  adding  a  deterministic  wide¬ 
band  signal  to  a  noise  sequence  and  then  estimating  the  signal  content 
at  the  loop  output  through  correlation  techniques,  such  as  Fourier  anal¬ 
ysis. 

If  one  accepts  the  model  assumption  that  the  combined  effects  of 
platform  movement  and  oscillator  dynamics  can  be  described  as  a  ran¬ 
dom  process  (paralleling  the  classical  requirement  that  the  PLL  have 
a  nonzero  signal  bandwidth)  then  our  approach  provides  an  objective 
method  of  comparing  the  performance  of  various  structures,  specifying 
an  optimum  performance,  and  testing  these  theoretical  results. 

4.  3  Time -Invariant  Nonlinear  Estimator 

4.  3. 1  Definition  of  Nonlinear  Filter  Class 

In  Figure  4-1  we  indicated  the  general  type  of  time -invariant  non¬ 
linear  filter  that  we  shall  consider;  the  use  of  a  Zero  Memory  Non- 
Linearity  (ZNL)  in  the  path  of  the  loop  error  signal  followed  by  a  time- 
invariant  linear  filter  to  generate  the  estimate.  The  selection  of  this 
filter  class  can  be  qualitatively  justified  by  considering  the  basic  signal 
processing  problem.  The  phase  message  is  a  slowly  varying  process, 
as  determined  by  platform  and  oscillator  dynamics,  whereas  the  noise 
is  uncorrelated  over  milliseconds  and  statistically  independent  over 

seconds.  Thus  we  expect  that  any  filter  matched  to  the  signal  will  have 

\ 

a  similarly  long  time  constant  and  it  would  be  desirable  to  limit  in  some 
way  the  large  amplitude  excursions  of  the  noise  prior  to  being  weighted 
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by  the  long  filter  time  constant.  Since  by  assumption  s?(k+l  |k)  is  small, 

the  loop  error  is  principally  composed  of  the  noise  component  and  it 

would  seem  logical  to  perform  the  nonlinear  amplitude  limiting  at  that 

point,  which  is  the  ZNL  location  of  Figure  4-1.  A  ZNL  followed  by  a 

linear  filter  with  memory  is  Zadeh's  class  (t^j)  nonlinear  filter  and  is 
52 

often  used  as  a  first  approximation  to  a  general  nonlinear  filter.  In 
section  4. 3.  5  we  shall  give  a  qualitative  argument  that  the  limiting  form 
of  the  optimum  ZNL  in  our  filter  provides  an  optimum  estimate  in  the 
maximum  likelihood  sense  and  our  experimental  results  will  bear  this 
out.  However,  in  the  following  sections  we  shall  be  concerned  with  anal¬ 
ysis  and  optimization  of  this  filter  class  which  we  have  specified  on  an 
intuitive  basis. 

4.  3.  2  Linearized  Analysis 

To  analyze  the  approximate  ZNL  class  filter  performance,  we  replace 
the  ZNL,  defined  by  f(z),  by  a  gain  normalized  ZNL,  f(z)/c,  where  c  is 
chosen  as 

s(k+l  Ik)  f(z(k+l)) 

c  -  - -  ;  - .  (4.17) 

s(k+l|k)2 

This  choice  allows  us  to  define  a  new  noise  component  at  the  output  of 
f(z(k+l)),  as 

n(k+l)  *  —  f(z(k+I))  -  s(k+l  [k)  (4. 18) 

c  ' 

which  is  orthogonal  to  the  signal,  s(k+l  Jk), 
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~n - 7  S  f(z)  -2 

s'—  f(z)  -  g  - - s  =  0. 

c 

s  f(z) 


The  gain,  c,  is  Ilooton's  equivalent  gain  and  the  retention  of  the  dis- 

4 

tortion  as  a  noise  term  is  an  extension  of  Booton's  method  due  to  Caron. 

This  type  of  linearized  analysis  of  a  given  nonlinear  control  system 

41  54 

appears  to  be  the  most  practical  technique  available  *  and  has  been 

used  by  several  workers  in  analyzing  PLL's  operating  in  gaussian 

45,  53 
noise.  * 

Recalling  our  assumption  that  s(k+l  jk)  is  small  relative  to  the  noise 
power,  R,  a  further  assumption  that  we  will  make  is  that  the  noise  se¬ 
quence  at  the  output  of  the  ZNL  is  uncorrelated,  n(k)  n(i)  =  0,  k*i. 

This  assumption  requires  that  the  uncorrelated  atmospheric  noise  se¬ 
quence,  y(k),  dominates  the  error  sequence,  sKk),  and  that  the  resulting 
ZNL  output  noise  sequence  n(k)  does  not  couple  significantly  back  to  the 
ZNL  input.  This  assumption  is  equivalent  to  the  qualitative  conditions 
listed  by  Smith  as  to  when  linearization  of  the  nonlinear  element  can 
be  used,  such  as  requiring  the  bandwidth  of  the  ZNL  input  to  be  signifi¬ 
cantly  larger  than  the  following  linear  filter.  If  this  assumption  is  valid, 
then  we  can  optimize  any  ZNL  type  filter  by  using  the  K-B  equations 
with  a  new  noise  power  n(k)  replacing  the  input  noise  power  R.  Fur¬ 
ther,  we  can  specify  the  optimum  filter  of  this  class  by  choosing  the 
ZNL  to  minimize  the  noise  power  n(k)  .  The  filter,  optimized  for  a 
given  ZNL,  can  then  be  used  to  check  the  n(k)  uncorrelated  assumption 
by  comparing  the  ratio  of  resulting  steady-state  estimation  error  (the 
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signal  input  to  the  ZNL,  s)  with  the  atmospheric  noise  input  to  the  ZNL. 
Our  results,  both  theoretical  and  those  derived  from  tests, indicate  that 
the  n(k)  uncorrelated  assumption  is  acceptable  for  most  problems  of 
practical  interest. 

If  we  define  a  signal-to-noise  ratio  improvement  factor 
I.  F.  *  SNRZnl  out/SNRZNL  IN' 

then  using  the  definition  of  a  gain  normalized  ZNL  and  (4. 17)  and  (4. 18) 
we  have 


R 

• 

r  2 

sf(z) 

~2 

-  —  2 

s 

s  •  1  (z)  -  s  f(z) 

(4.  19) 


If  the  assumption  that  n(kH)  is  an  uncorrelated  time  series  is  ralid,  then 
the  relations  between  the  I.  F.  and  the  performance  of  a  system  of  our 
class,  optimized  in  a  specified  ZNL,  is  given  functionally  as 


I.F.  =  g(R,s2  »  £^(00),  Pz(z)) 


/ 


r 

/ 


(4.  20) 
{4. 21) 


Sn(oo)  =  h(R/I.  F. ,  Signal  Model), 

/where  (4.  20)  depends  upon  the  first-order  pdf  of  z(k+l)  and  2^(00)  is  the 
K-B  steady-state  solution  (4. 15)  using  a  new  effective  noise  power, 
R/I.F.  Thus,  to  determine  the  performance  of  the  given  estimator 
class  we  can  study  the  I.F.  (4. 19)  for  various  types  of  ZNL's  and  predict 
the  resulting  performance  with  a  simultaneous  solution  of  (4.  20)  and 


(4.21). 
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If  we  assume  the  signal  generator  model  to  be  controllable  and  observ¬ 
able,  then  the  steady- state  solution  (4.21)  will  exist  and  it  will  be  a  mono- 

36 

tone  function  of  the  effective  noise  power  R/I.F.  The  minimum  mean- 
square  error  of  the  nonlinear  class  we  have  specified  (ZXoo)^.^)  will  then 
be  provided  by  the  ZNL  which  maximizes  the  Improvement  Factor. 

Hause  has  shown  that  the  I.  F.  is  maximized, at  a  given  R/  s  ,  by  the 
conditional  expectation  ZNL  which  provides  the  MMSE  estimate  of  s(k+l  |k) 
given  z(k+l)  and  is  defined  by 

n+  oo 

f(z)  =  \  T]  p  (x\,  z)  drj. 

J-oo  s  |  z 

.  ^2 

We  can  study  the  I.  F.  of  various  ZNL' s  versus  R/ s  and  compare 
them  to  the  optimum  MMSE  ZNL  as  a  measure  of  relative  ZNL  quality. 

If  the  n(k+i)  uncorrelated  assumption  is  valid,  we  can  also  determine 
the  simultaneous  solution  of  (4.  20)  arid  (4.  2  l)  for  a  specific  signal  gener¬ 
ator  model  from  the  evaluation  of  I.  F.  The  resulting  operating  point 
solution  will  then  provide  an  indication  of  the  validity  of  the  assumption 
in  terms  of  the  ratio  R/  s  . 


4.3.3  ZNL  Improvement  Factor 

To  study  the  ZNL  I.F.  for  a  given  R/s^  ratio,  we  must  evaluate 
(4. 19)  which,  in  terms  of  the  first-order  pdf's  is  given  by 

i2 

§ 


R 

l  ~2  J 


(E[E[s  f(z) |  sJJ  Y 

Li  mt 


s2  E[E[f2(z)|  s]  }  -  (Numerator) 

z  s 


I.F.  = 
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l/_p  (§)  J,  i P-JttS)  f(ii)  dtidil2 

s  g  z  y 

_  -  * 

s  2  /  p*(£)  /  p  (*]-£)  dqd£  -  [Numerator] 

s  5  y 

where  E[  •  ]  is  a  conditional  expectation  wrt,  x.  Interchanging  the  order 
of  integration  yields 

[  J  ^  S  PJh“£)  P J&  d^dTj]2 

_ s  y _ s _ 

s  2  /  f2(ri)  /  p  ($)  p  (r]-£)  d^dt]  -  (Numerator) 
s  s  y 


In  this  expression  p  ( • )  is  the  atmospheric  noise  pdf  and  p^(  • )  is  the 

^  s 

estimation  error  pdf.  To  proceed  further  we  must  know  the  latter  pdf 


which  is  formally  related  to  the  estimator  structure  and  non-gaussian 
observation  noise  by  the  Chapman-Kolmogorov  functional  equations.  We 


shall  assume  that  this  pdf  can  be  approximated  by  a  zero-mean  gaussian 
distribution,  an  assumption  that  is  consistent  with  previous  assumptions 
concerning  the  fact  that  s(k+l  |k)  is  small,  the  estimator  time  constant 
is  long  and  that  the  ZNL  will  attenuate  the  large  amplitude  excursions 
of  the  atmospheric  noise  so  that  no  single  noise  sample  win  dominate 
the  weighted  sum  within  the  time  constant  window. 

We  now  observe  that  the  interior  integration  in  the  left  term  of  the 
denominator  is  a  convolution  that  defines  p  ( • )  and  that  the  interior  inte- 
gration  of  the  numerator  is  proportional  to  the  derivative  of  a  gaussian 
pdf  convolved  with  the  noise  pdf.  Interchanging  these  latter  two  linear 


operations  yields 


Thus  the  right-hand  term  of  I.  F. ,  the  ZNL  SNR,  becomes 


SNR 


ZNL  “ 


[  /z  f<*l>  dTi]2 


T=r  4  Pz(t1>  dr]  ~  *•  4  PZW  dt1^2 

m2  “ 

s 


(4.22) 


where  p'  (ti)  represents  the  derivative  of  p  ( * )  wrt  its  argument, 
z  z 

The  ZNL  SNR,  (4.22),  is  maximized  by  the  ZNL  which  generates 
the  conditional  expectation  of  s,  which 

f(z)MMSE  =  J-  *>  P  ,  ^ z)  d1 
s  S  I  z 

4  ’i  Pv<z_ri)  P  Jn> 

s  y  s 


p  (r\)  as  above,  this  reduces  to 
s 

p  (z)  =  p  (£)  ®  p  (£).  (4.23) 

z  y  ^ 

J  s 

It  can  be  shown  using  calculus  of  variations  that  (4.  23)  does  maximize 
(4  22). 

Recalling  the  definition  of  the  atmospheric  noise  pdf  (3.  5)  we  see 
that  p,J,  • )  is  given  by 


Pz(z) 


Using  the  gaussian  assumption  for 


— o  P'(z) 
f  z  MMSE  -s  7^ 


IS 


■i 


T1  P  Pjn) 


z  s 


pz(z) 


dri 


dT! 


-131- 


f 


V 


pz(d  =  N(0,s2)  «  [q,N(o,^G)  +  P.n(o,4)  ®  p^lj 


“  q.N(°-'nG  +  S2)  +  P.N(°-'rnG+s2)  ®  PyH^’  (4‘24> 


For  small  values  of  z,  the  low  amplitude  gaussian  component  of  (4.  24) 
will  dominate  the  pdf  and  ^ Z^MIVISE  ^rom  ^4*  ^3)  will  be  a  s*raight  line 
passing  through  the  origin,  the  optimum  ZNL  estimator  for  gaussian 
signal  and  noise.  For  large  amplitude  values  of  z,  the  hyperbolic  com¬ 
ponent  of  Pyjj(z)  wiH  dominate  and  (4.  23)  will  yield  a  function  behavior 


Z2  Pk(z> 


-s 


a 


PZ<Z> 


_d _ l_ 

dz  m 
z 


(4.25) 


Thus  the  MMSE  ZNL  has  a  small  signal  linear  portion,  reaches 
a  point  of  maximum  value  and  decays  hyperbolically  for  large  z,  irre¬ 
spective  of  the  integer  m  value  of  the  atmospheric  noise  tail.  We  have 
evaluated  (4.24)  numerically  and  find  that  this  general  behavior  is  fol¬ 
lowed,  although  the  nature  and  position  of  the  peak  value  is  strongly 
dependent  upon  the  atmospheric  noise  weather  conditions  and  resulting 
pdf.  We  shall  not  consider  as  a  v^a^e  candidate  ZNL  to  use 

in  a  practical  estimator  since  one  has  to  know  a  great  deal  about  the 
specific  noise  pdf.  We  shall  use  the  SNRj^^  produced  by 
as  a  comparison  for  practical  ZNL's  to  be  discussed  below. 

4.  3.4  Improvement  Factor  of  Four  Types  of  ZNL's  in  Atmospheric 
Noise 

A  practical  ZNL  should  have  two  attributes,  achieve  near  optimum 
performance,  and  have  parameters  that  are  not  critically  sensitive  to 


i  i 
:  i 
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specific  noise  characteristics,  or  that  can  be  adjusted  in  a  noncritical 
manner.  The  first  and  most  obvious  candidate  ZNL  is  a  hard  limiter, 
defined  by 


f(z)  = 


1 

-1 


z  5*  0 
z  <  0. 


(4.2  ) 


The  hard  limiter  in  this  position  in  a  sampled  data  PLL  can  be  simply 
realized  by  hard  limiting  the  received  rf  signal  at  the  output  of  the  antenna 
bandlimiting  filter  and  sampling  the  resulting  two-state  waveform  at  the 
time  of  the  reference  carrier  (VCO  output)  zero  crossing.  This  is  not 
exactly  the  same  type  of  system  that  one  obtains  by  hard  limiting  both 
received  signal  and  carrier  (often  termed  bandpass  limiting),  multiply- 

45 

ing  and  then  processing  the  resulting  waveform  with  a  lowpass  loop  filter. 
One  continuous  time  dual  of  the  hard  limiter  that  we  are  considering  is 
the  use  of  the  limiter,  after  lowpass  filtering  to  recover  the  baseband 
phase  error  signal,  i.  e. ,  a  bang-bang  control  system  application.  The 
hard  limiter  ZNL  does  not  have  a  parameter  explicitly  associated  with 
its  transfer  function,  however,  one  must  know  the  value  of  c,  (4. 17), 
the  signal  suppression  of  the  limiter, so  that  the  optimum  loop  gain, 

K/c,  can  be  specified  for  optimum  performance  for  a  given  signal  model 
and  noise  condition. 

Two  other  types  of  ZNL's  that  have  been  proposed  in  the  literature 
are  clipping  and  hole -punching,  defined  by 


z,  |z|  «Zt 

zr  Z<zt 


f(z) 


hP 


z  ^Z, 


0,  z  >Zj.. 


(4.27) 
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These  two  ZNL's  are  shown  in  Figure  4-2,  along  with  the  limiter  and 
the  qualitative  behavior  of  the  MMSE  ZNL  as  numerically  evaluated 
from  (4.  23).  The  clipping  and  hole-punching  ZNL's  have  a  single  thresh¬ 
old  parameter  to  define  their  transfer  function,  and  we  must  either  carry 
this  as  an  additional  parameter  in  the  analysis  or  specify  them  as  a 

function  of  the  noise  characteristics.  We  have  adopted  the  latter  and 

42 

used  a  suggestion  by  Feldman  that  the  threshold  be  set  so  that  a  fixed 
percentage  of  the  samples,  z(k+l),  fall  in  the  nonlinear  region.  This 
behavior  can  be  easily  controlled  by  a  long-time  constant  digital  servo 
loop  that  simply  counts  samples  exceeding  the  threshold  and  then  adjusts 
the  threshold  in  the  correct  direction  at  the  end  of  each  averaging  period. 
Thus  the  technique  is  not  dependent  upon  a  specific  knowledge  of  noise 
conditions;  however,  we  must  evaluate  the  I.  F.  for  various  noise  con¬ 
ditions  to  determine  if  the  performance  using  this  technique  is  insensi¬ 
tive  to  noise  conditions. 

Substitution  of  (4.  23)  into  (4.  22)  yields  the  maximum  I.  F.  for  a 
given  noise  pdf. 


(I.  F .) 


max 


(4.  28) 


Substitution  of  the  hard  limiter  characteristic  results  in  a  simple  expres¬ 
sion  for  the  I.  F. 


I.F. 


h.  1. 


1 


«<2  i  2 /  _ \ 

s  4p  (0) 


(4.29) 
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Figure  4-2.  Various  ZNL  Transfer  Functions 
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These  I.  F.'s,  including  the  clipper  and  hole-puncher,  were  numer¬ 
ically  evaluated  for  all  four  atmospheric  noise  conditions  of  Chapters  2 

and  3,  and  for  gaussian  noise.  The  results  are  shown  in  Fig.  4-3  through4-5 

$ 

representative  percentages  and  the  hard  limiting  and  optimum  cases. 

Inspection  of  these  results  suggests  the  following  comments: 

a)  The  ZNL  I.  F.  plots  show  only  a  15%  clip  percentage.  Perform¬ 
ance  was  evaluated  at  5%  increments  from  5%  to  50%  and  this 
percentage  was  selected  as  providing  the  best  I.F.  over  the 
widest  possible  range  of  conditions.  We  note  that  the  15%  clipper 
provides  near  optimum  performance  over  the  widest  range  of 
R/s  value  of  any  practical  ZNL.  It  is  thus  insensitive  to  both 
the  range  of  input  signal -to -noise  as  well  as  being  relatively  in¬ 
sensitive  to  the  type  of  atmospheric  noise. 

b)  Hole  punching  is  clearly  more  sensitive  than  clipping  to  both  the 
type  of  atmospheric  noise  and  the  R/s  ratio.  While  it  is  pos¬ 
sible  to  select  a  good  hole  punch  percentage  for  a  specified  oper¬ 
ating  point,  significant  changes  in  conditions  will  render  this 

selection  a  poor  one.  Similar  results  for  hole  punching  have 

1 Q  51 

been  reported  by  Kapp  and  Griffiths  for  the  signal  detection 
problem.  They  plotted  a  measure  of  detection  performance 
against  a  single  noise  parameter  and  found  hole  punching  exhib¬ 
ited  a  much  sharper  performance  peak  than  clipping  (they  did  not 
use  our  percentage -threshold  technique). 

.Z 2 

c)  Hard  limiting  performs  quite  well  for  most  regions  of  low  R/  s 

*Notes:  The  straight  lines  on  these  plots  are  the  functions  (4.21)  for  two 
signal  models  to  be  given  in  section  4.  5;  *\/r/  s  2  »  [Z^/R]  1/2. 
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Figure  4-3.  ZNL  Improvement  Factor  for  Frontal  and  Tropical  Noise 
Conditions 
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Figure  4-4.  ZNL  Improvement  Factor  for  Quiet  and  Quiet-Night 
Noise  Conditions 
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Figure  4-5.  ZNL  Improvement  Factor  for  Gaussian  Noise  and 
Equivalent  Linear  Gain 
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ratio  where  we  expect  operating  solutions  of  the  estimator  struc¬ 
ture  to  validate  the  n(k)  uncorrelated  assumption.  In  regions 
where  the  hard  limiter  performance  decreases  rapidly,  this  un¬ 
correlated  assumption  will  also  break  down  since  the  correlated 
signal,  s,  will  cause  a  significant  portion  of  the  hard  limiter 
distortion  (noise  n(k)).  This  will  in  general  make  the  actual  per¬ 
formance  less  than  predicted  by  the  combination  of  the  I.F. 
curves  and  K-B  equations  which  assume  that  n(k)  is  uncorrelated, 
d)  Figure  4-5  also  shows  the  variation  of  c  (4.  17),  with  noise  con¬ 
dition,  for  hard  limiting  and  1 5%  clipping.  Following  the  develop¬ 
ment  of  4.  3.  3  it  can  be  shown  that  this  parameter  is  given  by 

c  =  -Jz  f(q)  p'z(ri)  dq,  ch  1  =  2P2^ 

which,  as  noted,  is  Booton's  equivalent  gain. 

Since  the  hole  punching  ZNL  has  distinct  performance  disadvantages 
compared  to  the  clipper  and  is  only  marginally  easier  to  implement,  we 
have  not  considered  it  further  in  this  research.  In  sections  4.  5  and  4.  8 
the  ZNL  performance  analysis  developed  here  will  be  combined  with  a 
scalar  and  two-state  signal  model  to  analyze  ZNL  estimator  performance. 
The  predicted  performance,  based  on  the  linearized  analysis  given  here 
will  be  found  to  agree  well  v/ith  tests  using  noise  sample  records  and  a 
simulated  signal  process. 

4.3.5  Optimality  of  ZNL  Estimator  Class 

The  ZNL  I.F.  curves  given  above  share  a  common  characteristic 

~  O 

in  that  for  decreasing  s  /R  ratio  they  approach  a  constant  value  for  any 


ZNL  characteristic.  This  value  is  determined  by  _Lim 


For  the  MMSE  ZNL  we  see  that  this  results  in 
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pz(z)  =  P  (z>. 


f(z) 


MMSE 


-s2 


PL(z> 

Pz<z> 


_  P'Jrt  - 

~ >2  "  ~  2  d  ,  t-  /  %  i 

-s  —  =  's  dilntPy(z>1' 

py  2 


(4.31) 


Several  workers^*  ^  have  shown  that  this  ZNL  (4.  31),  followed  by 

a  matched  filter,  is  the  asumptotically  optimum,  weak  signal  detector, 

for  known  signals  in  non-gaussian  noise  where  noise  samples  are  statis- 

37 

tically  independent.  It  can  be  shown,  combining  Van  Trees'  develop¬ 
ment  of  a  maximum  likelihood  estimate  of  a  discrete  parameter,  s,  with 
44 

Antonov's  technique  for  a  Taylor  series  expansion  of  the  likelihood 
ratio,  that  an  asymptotic  optimum  estimate  of  a  parameter  s  is  given 
by  the  ZNL  of  (4.  31)  followed  by  an  averaging  function.  Furthermore, 
the  large  amplitude  form  of  the  noise  pdf  (3.  14a)  satisfies  Capon's  con¬ 
dition^ 


0  <  f-x  (PyW2/PyW  dll  <  00 * 

A 

to  ensure  that  the  resulting  estimate  s  tends  to  gaussianity.  We  can 
argue  that  our  assumed  ZNL  estimator  structure,  feeding  back  the  esti¬ 
mate  of  s(k+l  Jk)  to  maintain  s  small,  should  provide  an  optimum  esti¬ 
mator  structure  for  a  long  observation  time  where  we  can  "open"  the 
feedback  loop  and  consider  the  estimation  of  s(k+l)  to  be  a  sequence  of 
discrete  parameter  estimation  problems.  We  have  not  obtained  a  formal 
proof  that  the  optimum  ZNL  class  filter  is  aymptotically  equivalent  to 
an  unstructured  optimum  due  to  the  difficulty  of  including  a  model  for 
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the  signal  dynamics  within  the  loop.  However,  we  shall  find  that  the  near 
optimum  estimate  obtained  with  the  use  of  the  pilot  channel  information 
and  the  time -varying  gaussian  noise  model  never  improves  upon  the 
error  performance  predicted  for  MMSE  ZNL.  We  suspect  that  the 
MMSE  ZNL  estimator  provides  the  optimum  estimate  of  all  filters  when 
sufficiently  dense  sampling  is  employed  and  the  noise  samples  are  statis¬ 
tically  independent.  For  the  case  of  non-independent  samples,  as  is 
found  in  several  atmospheric  noise  conditions,  the  optimum  performance 
will  be  less  than  the  MMSE  ZNL  estimator. 

4.4  Linear  Time-Varying  Estimator 
4.4. 1  Use  of  Joint  Channel  Noise  Model 

We  demonstrated  in  Chapters  2  and  3  that  a  useful  model  for  a  band- 

limited  atmospheric  noise  waveform  is  a  time-varying  gaussian  process. 
3fi  37 

It  is  well  known  ’  that  the  general  time-varying  Kaman-Bucy  filter 
provides  the  MMSE  estimate  of  a  rational  spectrum  signal  process  addi- 
tively  combined  with  a  time -varying  gaussian  noise  process  where 

E[y(t)  y(t+T>]  =  R(t)  6(t-r). 

The  experimental  sample  records  on  which  we  have  based  our  noise 
model  utilized  the  pilot  channel  envelope  as  a  measure  of  the  rf  noise 
variance.  Equation  (3.  23)  gave  the  conditional  estimate  of  the  noise 
power,  conditioned  by  the  pilot  channel  envelope  value  v  ,  as 

2 

2  -  1  nH 

y|vo,x=l  "  m  2 

°vH 


2  2 

v  +m<r  „ 
o  vH 


(3.23) 

repeated 
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for  noise  excursions  above  the  background  gaussian  noise  process.  To 
utilize  the  pilot  channel  envelope  as  a  measure  of  the  time -varying  rf 
noise  power  we  shall  implement  the  following  algorithm: 


R(k)  =  <r 


^k 


nGJ 


1 

m 


nR 

2 

rvH 


[4 


<  v. 


>  V, 


(4.32) 


The  threshold  value  of  the  envelope,  vt,  was  selected  as  the  intersection 

of  the  relation  (3.  23)  with  the  it  q  value  as  seen  on  the  model-data  com- 

2  2  2,2 

parison  Figures  3-1,  3-2.  Therefore,  vt  =  monGtrvH<  °nH* 

Both  the  joint  channel  model  and  the  experimental  observations  show 
that  the  pdf  of  the  rf  noise  samples,  y(k),  is  not  gaussian  when  condi¬ 
tioned  on  the  pilot  channel  envelope  even  though  the  model  utilizes  a 
noise  process  in  each  channel  which  is  gaussian  when  conditioned  on  the 
frequency  independent  random  variable  "a".  If  the  model  were  an  exact 

description  of  the  joint  channel  noise  relations,  one  could  use  a  better 
2 

estimator  of  c  i  than  provided  by  the  pilot  channel  envelope.  However, 

y  I  vo 

the  envelope  was  used,  based  on  empirical  evidence  that  it  provided  the 
highest  correlation  coefficient  with  |  R,I:‘  |  ^  and  hence  we  must  adapt  our 
use  of  the  joint  channel  model  and  time -varying  filter  to  our  experimental 


records.  To  consider  the  effect  of  the  non-gaussian  pdf  of  y|vQ  we  eval¬ 
uated  the  MMSE  ZNL  Improvement  Factor,  as  described  in  section  4.  3, 
for  the  conditional  distribution  (3. 11)  with  m  =  1  and  found  that  a  1. 35  dB 
improvement  was  possible  in  estimating  s  given  z.  To  reflect  this  pos¬ 
sible  improvement  in  an  ad  hoc  manner  in  the  time -varying  estimator. 
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we  utilized  a  clipper  ZNL  with  threshold  set,  at  each  iteration,  to  plus 
and  minus  three  sigma  of  the  sum  of  ihe  estimate  error  and  noise  for 
that  iteration.  This  algorithm  will  be  given  in  sections  4.  5  and  4. 6 
when  we  consider  the  actual  simulations.  Thus  the  algorithm  that  we 
shall  use  to  specify  the  time-varying  power  level  of  the  atmospheric 
noise  sample  sequence  will  not  result  in  a  precisely  MMSE  estimate 
the  signal  process  since  the  conditional  noise  distribution  is  not  gaussian. 
The  above  computation  and  the  ad  hoc  correction  scheme  should  result 
in  performance  that  is  within  1  dB  of  optimum. 

4.  4.  2  Performance  of  the  Time-Varying  Linear  Estimator 

The  qualitative  effect  of  the  time-varying  power  level  on  the  estima¬ 
tor  can  be  understood  by  referring  to  the  K-B  sampled  data  equations 
(4. 12),  and  (4. 13)  and  (4. 14).  The  time-varying  gain,  K(k+1),  responds 
to  each  R(k+1)  value  in  an  inverse  manner,  thus  attenuating  samples, 
z(k+l),  that  are  caused  by  large  noise  excursions  beyond  the  background 
noise  level.  To  the  extent  that  we  can  associate  an  expected  magnitude 
with  the  R(k+1)  estimate,  |z(k+l)|  =  a*  R(k+1)1/^,  we  see  that  the  R  *(k+l) 
portion  of  (4. 12)  provides  a  transfer  relation  quite  similar  to  the  MMSE 
ZNL  (4.  25);  it  is  linear  for  R(k+1)  <  vt  and  decays  hyperbolically  beyond 
this  point. 

The  second  time-varying  term  of  (4.  12)  introduces  a  dynamic  be¬ 
havior  into  the  attenuator.  This  term,  2(k+l  |k+l)  reflects  the  growth 
in  uncertainty  of  the  signal  state  at  each  iteration  (4.  13)  and  the  decrease 
in  uncertainty  caused  by  optimal  weighted  observations  (4.  14).  When  a 
significant  number  of  sequential  observations  are  heavily  contaminated 
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by  noise,  the  growth  of  2(k+l  |k+l)  will  increase  the  weighting  of  subse¬ 
quent  samples  that  have  a  small  noise  contamination.  This  increase 

2 

will  be  above  the  quiescent  value  caused  by  R  .  (k+1)  =  <r  ~  and  will  be 

^  mm  nu 

36 

a  stochastic  manifestation  of  the  well-known  fact  that  a  K-B  filter  has 
a  step  response  that  is  "faster”  than  a  time -invariant  filter  equal  to  the 
steady- state  K-B  filter. 

Applying  this  description  to  our  atmospheric  noise  model  we  see 
that  for  quiet  conditions  where  there  is  no  significant  multiple  discharge 
activity  and  the  noise  samples  are  statistically  independent  beyond  sev¬ 
eral  milliseconds,  the  time -varying  estimator  should  have  a  performance 
quite  similar  to  the  ZNL  estimators.  For  heavier  noise  conditions  such 
as  tropical  and  frontal,  where  the  samples  are  not  statistically  indepen¬ 
dent  due  to  the  multiple  discharge  activity,  we  would  expect  that  the 
time-varying  filter’s  optimal  compensation  for  this  via  2(k+l  jk+1)  should 
result  in  improved  performance  over  the  ZNL  estimator. 

To  predict  the  performance  of  the  time-varying  filter  we  must  de¬ 
termine  the  average  value  of  2(k+l  |k+l).  From  (4. 13)  and  (4. 14)  we 
see  that  this  requires  that  we  determine  the  first-order  output  statistics 
of  a  nonlinear  difference  equation  driven  by  a  non-gaussian,  colored 
process,  R(k),  as  defined  by  (4.  31).  There  are  no  practical  solution 
techniques  for  this  problem  short  of  linearization  of  the  equation  and 
approximation  of  the  driving  process  with  a  gaussian  process.  Since 
the  actual  driving  process  is  very  non-gaussian,  this  latter  step  would 
seem  to  make  the  approximation  of  questionable  value.  The  one  prac¬ 
tical  answer  that  can  be  obtained  is  to  place  an  optimistic  lower  bound 
on  the  error  variance  by  considering  the  steady-state  error  for  a  K-B 
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filter  with  constant  noise  power  Rmin  =  <rnQ>  i.e. ,  ignoring  the  effect 
of  all  the  non-gaussian  excursions.  Thus 

~T.  V.(c°)  ^“RMIN(co)* 

In  the  two  simulation  examples  that  we  have  explored,  a  tighter  bound  is 
found  to  be  given  by  the  predicted  MMSE  ZNL  estimator  error.  If  this 
ZNL  filter  does  provide  an  asymptotic  optimum  estimator  for  statisti¬ 
cally  independent  noise  and  long  observation  times  (dense  sampling  wrt 
filter  time  constant)  as  we  suspect,  then  it  follows  that  shorter  observa¬ 
tion  intervals  and  statistically  dependent  noise,  though  processed  in  an 
optimal  manner  by  a  time-varying  filter,  cannot  result  in  better  perfor¬ 
mance  than  the  asymptotic  optimum. 

4.  5  Test  of  Estimators  for  Scalar  Signal  Model 
4.  5. 1  Operating  Equations 

The  first  signal  model  we  shall  consider  is  a  scalar  case  with  F  =  0, 
which  is  a  single  integrator.  This  model  is  useful  for  certain  navigation 
applications  where  the  mobile  platform  uses  a  precision  reference  oscil¬ 
lator  and  the  single  integration  of  the  white  guassian  driving  process 

represents  phase  instabilities  in  the  transmitting  and  receiving  oscil- 
38 

lator  and  a  random  velocity  description  of  the  mobile  platform  move¬ 
ments.  The  sampled  data  description  of  the  signal  generator  is  seen 
to  be  (equations  (4.  l)-(4.4)) 

$(AT,0)  =  1,  E[u(k)  u(k)T]  =  QAT. 

The  sampled  data  estimator  equations  can  be  solved  directly  for  steady- 


state  to  yield 
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2n(oo)  =  - 


QAT  +  V  (QAT r  +  4EQAT 


«  n/RQAT, 


R  »  QAT 


(4.33) 


QAT 


S(k+l|k+l)=  1 


QAT 


s(k|k)  + 


'QAT 


r(k+l),  R  »  QAT. 


(4.34) 


The  estimator  is  thus  a  filter  with  a  sampled  exponential  impulse  response. 
To  relate  our  model  to  classical  design  techniques  (section  4.  2.  3)  one 
would  specify  a  noise  power  level  and  filter  time  constant,  AT,  and 
determine  the  required  signal  model  driving  process  variance,  QAT, 
from 


=  .37, 


1/Nf  2 

QAT  =  R [l -.37  M  . 


To  determine  the  steady-state  estimate  error  of  the  ZNL  estimator 
we  utilize  (4.  32),  the  definition  of  the  Improvement  Factor  and  R  to  ob- 


(RQAT  _  RQAT  2SNR  •  QAT 

t.F.  •  '-F-=  =  (^TTr)4 


(4.35) 


where  we  have  used  £  j  j  «  s  ^(k+1  |k).  The  equation  (4.  35)  has  been 
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plotted  on  the  ZNL  I.F.  curves,  figures  4-3  through  4-5  for  represen¬ 
tative  values  of  the  SNR  •  Q  product,  with  AT  =  1  for  convenience.  The 
intersection  of  these  lines  with  the  ZNL  I.F.  curves  then  represent  a 
graphical  solution  of  the  two  conceptual  relations  (4.  20)  and  (4.  21).  For 
a  given  set  of  model  parameters,  Q,  and  SNR,  one  determines  the  esti¬ 
mator  operating  point  'sT2~j  /r,  at  this  intersection,  which  can  be  con¬ 
verted  to  estimation  error  variance  in  radians  using  R  =  l/2  SNR.  The 
operating  point  for  the  linear  time -invariant  system  is  the  intersection 
of  the  Q  •  SNR  line  with  the  horizontal,  0  dB  I.F.  axis. 

4.  5.  2  Simulation  Technique 

Figure  4-1  shows  the  general  estimator  block  diagram  plus  the  inter¬ 
action  with  the  atmospheric  noise  sample  record.  This  sample  record 
is  described  in  Chapter  2  and  consisted  basically  of  10-second  bursts  of 
nearly  simultaneous  samples  of  the  rf  noise  in  frequency  channel  1  and 
the  envelope  in  channel  2.  The  sampling  rate  was  1  kHz  and  with  total 

5 

length  of  approximately  10  samples.  Only  the  magnitude  of  the  rf  noise 
samples  was  recorded  and  this  sample  record  was  converted  to  a  bipolar 
record  using  a  random  number  generator  to  assign  a  sign  bit  to  each 
sample.  This  provided  an  uncorrelated  time  series  as  the  results  of 
Chanter  2  indicated  was  appropriate  for  sampling  intervals  greater  than 
the  bandlimiting  filter  time  constant.  A  computer  generated  time  series 
of  statistically  independent,  normal  samples  with  unit  variance,  was 
written  "alongside"  the  atmospheric  noise  record  to  provide  the  signal 
model  driving  process.  As  indicated  on  Figure  4-1  the  noise  and  signal 
samples  were  scaled  by  the  appropriate  factor  as  they  were  read  from 
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the  magnetic  tape  for  each  simulation  run  to  provide  the  specified  Q 
and  R  levels  for  that  simulation. 

For  the  linear  time -invariant  estimator,  the  loop  parameters  are 
given  by  (4.  33).  For  the  two  ZNL  types  that  were  tested,  15%  clipping 
and  hard  limiting,  (4.  34)  was  used  with  R  given  by 


2SNR  •  I.  F.  (Solution) 


where  I.  F.  (solution)  indicates  the  I.  F.  at  the  predicted  estimation  error 
operating  point  in  Figures  4-3, 4-4.  In  addition,  the  c  value  curve. 

Figure  4-5,  was  used  by  entering  at  the  predicted  \l  2},/R  operating  point 
and  reading  the  c/4r  value.  The  factor,  l/c,  was  applied  following  the 
ZNL  to  maintain  average  signal  gain  through  the  device  at  unity.  In  the 
case  of  1 5%  clipping,  the  threshold  yielding  1 5%  clipping,  which  was 
computed  in  the  ZNL  performance  analysis,  was  used  and  rescaled  to 
the  R  value  of  each  simulation.  This  threshold  can  be  easily  set  by  a 
long-time  digital  servo  loop  as  previously  noted.  This  approach  was  not 
used  here  to  allow  use  of  the  entire  sample  record  for  testing  of  the  per¬ 
formance  analysis.  To  recognize  the  fact  that  the  ZNL  estimator  analysis 
was  only  approximate,  three  simulations  were  run  in  parallel  for  the 
clipper  and  hard  limiter  with  slightly  smaller  and  larger  loop  bandwidths 
than  predicted  (realized  with  a  change  in  c  value)  to  verify  that  the  com¬ 
posite  gain  term,  yielding  minimum  estimation  error,  was  correctly 
given  by  the  analysis. 

2 

For  the  time-varying  filter,  two  parameters,  crnQ  A  Rmin  and  vt, 
were  entered  for  each  noise  record  to  describe  the  experimental  curve. 


2 

j  v  ^).  These  were  related  to  R(k+1)  by  the  algorithm  of  section4.4.1. 
The  time-varying  value  of  K(k+1)  wus  then  computed  from  (4.  12)-(4. 14) 
with  2(0  1 0)  =  0  although  this  initial  condition  decayed  so  quickly  with 
respect  to  the  sample  record  length  as  to  have  no  effect.  The  time- 
varying  clipper,  recognizing  the  fact  that  the  conditional  rf  noise  distri¬ 
bution  was  not  exactly  gaussian  (section  4.4. 1),  was  implemented  as 


Z(k+1)  = 


Z(k+1), 

.±ZT’ 


|  Z(k+1)  |  *5  3  y  2(k+l  |  k)  +  R(k+1)  &  Zt 
Z(k+i/  <  ±Z,p. 


(4.36) 


The  complete  simu~’tion  then  consisted  of  four  types  of  PLL  esti¬ 
mators  being  dri .  en  by  the-  same  sample  record  of  simulated  signal  and 
recorded  atmospheric  noise  samples.  The  resulting  estimation  error 
was  evaluated  at  each  iteration  as  seen  in  Figure  4-1,  with  this  error 
used  to  plot  the  time  evolution  of  the  error  and  determine  the  error  vari¬ 
ance.  The  results  of  the  scalar  simulation  for  all  noise  conditions  is 
given  below. 

4.  5.  3  Simulation  Results  for  Scalar  Model 

Four  plots  of  the  standard  deviation  of  the  estimate  error  are  shown 
in  Figure  4-6  for  the  scalar  model  simulation.  The  effective  noise  im¬ 
provement  scale  is  computed  from  (4.  35).  The  upper  right  plot  repre¬ 
sents  the  widest  effective  estimator  bandwidth  tested  and  we  see  that  the 
time-varying  filter  performs  slightly  better  under  all  conditions,  while 
the  hard  limit  ZNL  is  definitely  the  poorest  estimator  form.  Referring 
to  Figure  4-3,  the  hard  limit  I.F.  for  frontal  noise  conditions  falls  below 
that  of  clipping  at  this  Q  SNR  product  and  this  fact  is  reflected  in  the 
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Figure  4-6.  Standard  Deviation  of  Estimation  Error  for  Scalar  Signal 
Model 
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simulation  results.  The  other  noise  conditions,  where  the  hard  limiter 
is  poorer  at  this  QSNR  product,  predict  similar  large  separations  in 
the  performance  between  hard  limiting  and  clipping  and  this  is  found  in 
the  simulations. 

As  the  estimator  bandwidth  is  reduced  by  decreasing  Q  and  SNR 
(upper  left,  lower  left,  lower  right  in  that  order)  we  see  that  the  per¬ 
formance  of  hard  limiting  and  clipping  draws  closer  together,  both  pre¬ 
dicted  and  tested.  In  the  lower  left,  a  crossover  between  clipping  and 
hard  limiting  is  predicted  and  found  in  the  test  results.  This  is  clearly 
seen  in  Figure  4-3  for  frontal  noise  conditions.  The  time-varying  opti¬ 
mum  filter  follows  the  best  ZNL  performance  and  generally  exceeds  it 
■  slightly  as  the  MMSE  ZNL  performance  predicts.  The  least  separation 

between  the  best  ZNL  tested  and  time -varying  is  found  in  the  narrowest 
,  estimator  bandwidth  test,  lower  right.  Here,  the  effects  of  noise  power 

level  correlation  are  erased  by  the  long  processing  time  involved.  All 
that  is  required  is  the  proper  ZNL  noise  limiting  in  the  estimator,  the 
dynamic  performance  advantage  of  the  time-varying  filter  is  averaged 
out  by  the  processing  time. 

A  plot  of  the  time  evolution  of  the  estimation  error  is  shown  for  two 
conditions  in  Figure  4-7.  In  the  upper  traces,  the  time-varying  filter 
is  seen  to  be  clearly  superior  while  in  the  lower  plot  the  hard  limiter 
1  appears  somewhat  better.  Both  of  these  traces  are  plotted  on  the  same 

time  scale,  illustrating  the  differences  in  estimator  bandwidth.  We  shall 
,  defer  comments  as  to  the  implications  of  these  results  for  receiver  de¬ 

sign  until  section  4.  7.  The  results  do  appear  to  varify,  at  least  for  the 
scalar  signal  model,  the  applicability  of  the  linearized  ZNL  estimator 
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analysis  and  this  application  of  the  joint  channel  noise  model  with  attendant 
concept  of  a  time-varying  gaussian  noise  process. 


4. 6  Performance  of  Estimators  for  Two-State  Signal  Model 
4.6.1  Operating  Equations 

The  most  commonly  used  signal  model  for  the  navigation  problem 
is  two-state. 


a  double  integrator  model.  Physically,  this  form  of  signal  generator, 
models  platform  movement  as  a  random  acceleration  and  allows  for  both 
phase  (Sj(t))  and  frequency  (Sgft))  offsets  between  the  transmitter  and 
receiver  oscillators.  The  sampled  data  equations  ((4. 10)-(4. 14))  for  this 


(4.37) 

0 

QAT 


The  sampled  data  steady-state  covariance  equations  ((4. 13),  (4. 14)) 
lead  to  fourth-order  algebraic  equations  for  this  model  and  they  cannot 
be  solved  analytically.  We  therefore  solve  the  steady-state  continuous 
time  equations,  (4.9),  to  obtain 


The  steady-state  estimation  error  and  gain  coefficients  for  the 
equivalent  sampled  data  system  is  obtained  by  substituting  the  results 
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of  section  4.  2.  2  into  (4.  39).  For  the  phase  estimate  variance  we  obtain, 
with  the  introduction  of  the  ZNL  Improvement  Factor, 


Sn(oo)  =  nT2 


RAT 

I.F. 


3/4  1/4 

Q  *  > 


I.F.  = 


(8Q  •  SNR  •  AT  ) 


3.1/3 


N2u/R) 


8/3 


(4.41) 


Expression  (4.41)  is  plotted  on  the  ZNL  I.F.  curves.  Figures  4-3, 
4-5,  for  a  representative  range  of  Q  •  SNR  products,  with  AT  again  set 
equal  to  one  for  convenience.  These  plots  show  that  the  two-state  esti¬ 
mator  operating  point  lines  have  less  slope  than  in  the  scalar  case, 
which  is  due  to  the  coupling  of  errors  between  the  phase  and  frequency 
states  (Sj(k),  s2(k)).  We  also  note  that  a  large  spread  of  parameter 
values  is  more  concentrated  in  a  region  of  the  I.F.  curves  indicating 
that  there  should  be  less  change  in  relative  estimator  performance  for 
different  ZNL's  with  this  signal  model.  The  sampled  data  estimator 
for  the  time -invariant  simulations  is  given  by 


Sj(k+1  |k+l) 

1  AT 

^  /8Q  •  SNR  I.F. 
V  AT 

= 

I(kjk)  4- 

s2(k+l  |k+l) 

0  1 

/ 2Q  •  SNR  I.F. 

V  AT 

•  -£■  f(r(k+l)  -  Sj(k+1 1 k)); 


(4.42) 
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c  =  1,  I.F.  =1,  f{*)=  *  for  linear;  c,I,F.,f(  •)  —  ZNL  solutions  for  nonlinear. 


Using  the  two-state  sampled  data  signal  model,  (4.37)  and  (4.38) 
and  the  time -varying  covariance  matrix  recursions  (4. 13)  and  (4. 14) 
we  can  determine  the  sampled  data  covariance  recursion  as 


r(k+l|k)  = 


211  +  2ATE12  +  At2z:22 
S21  +  S22AT 


212  +  ATS22 
S12  +  QAT 


2..  =  2..(k), 


r(k+l|k)  =  2(k+l  k). 


rnR(k+i) 

rn  +  R(k+1) 


2(k+l  |  k+1)  = 


r21R(k+D 
ru  +  R(k+1) 


r12R(k+i) 

rn  +  R(k+1) 

rnr22  ~  ri2  *  r22R(k+I) 
rn  +  R(k+1) 


(4.43) 


(4.44) 


r..(k+l|k). 


=  r 


21* 


Note  that  here  we  are  using  the  actual  sampled  data  model  based  on  the 
physical  continuous  time  system  and  hence  use  QAT  and  R  as  the  signal 
driving  process  and  noise  variance,  respectively.  The  recursions  (4.43) 
and  (4.44),  which  cannot  be  solved  analytically  for  the  steady  state  2(oo) 
value,  are  required  in  the  time-varying  filter  to  compute  the  gain  K(k+1). 
This  gain  is  computed  for  each  sample  iteration  from  (4.43),  (4.44), 
and  (4. 12).  To  verify  that  the  approximate  sampled -data/continuous 
relations,  in  conjunction  with  (4.  39),  yield  the  value  to  which  (4.43) 
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and  (4. 44)  converge,  they  were  iterated  nr."ierically  using  a  sampling 
interval 


where  N  is  the  number  of  samples  per  unit  of  filter  time  constant.  The 
actual  iterated  values  of  the  2(oo)  matrix  agreed  to  within  1%  of  those 
given  by  substituting  R  =  RAT  into  (4.  39)  for  N  >  100. 

''T 

4,6.2  Simulation  Technique 

The  simulation  technique  used  for  the  two-state  signal  model  was 
essentially  that  given  for  the  scalar  case.  The  estimator  recursion  was 
given  by  (4.  42)  with  the  appropriate  ZNL  parameters  determined  by  the 
graphical  solutions  of  (4.41)  with  the  I.  F.  curves.  The  same  mapping 
from  E(k+1)  to  R(k+1)  was  used  for  the  time- varying  filter,  with  the 
estimator  equations  given  by  (4.43)  and  (4.44)  plus  the  general  K-B 
equations.  The  time-varying  clipper  was  also  used  as  per  (4.36), 

4.6.  3  Simulation  Results  for  Two-State  Model 

The  standard  deviation  of  the  phase  estimate  error  is  shown  in 
Figure  4-8,  in  the  same  format  as  used  for  the  scalar  model  results. 
The  three  plots  are  for  decreasing  estimator  bandwidth,  upper  left, 
lower  left,  and  right,  respectively.  In  the  first  plot,  we  note  the  mag¬ 
nitude  of  the  error  is  such  that  we  cannot  conclude  that  the  small 
s(k+l  |k)  assumption  holds  and  we  must  interpret  this  as  merely  an 
extension  of  the  linear  model  to  determine  the  effects  of  bandwidth  ver¬ 
sus  the  noise  conditions  and  estimator  structures. 
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A  general  comment  that  applies  to  all  of  these  results  is  that  the  two- 
state  model,  with  its  coupling  between  phase  and  frequency  state  errors, 
produces  many  instances  of  significant  deviations  between  linearized  pre¬ 
dictions  and  actual  results.  We  attribute  this  to  the  state  coupling  and 
the  effect  of  multiple  discharges  in  the  atmospheric  noise  record  that 
reduce  the  information  content  of  long  sample  sequences.  In  all  of  the 
frontal  noise  simulations  the  optimum  time-varying  filter  performs  at 
least  3  dB  better  due  to  its  improved  dynamic  response  and  optimal 
weighting  of  uncontaminated  samples.  This  improvement  is  also  seen 
in  the  tropical  noise  record  recorded  in  the  afternoon  (right-hand  set  of 
points)  as  compared  to  the  morning.  These  two  records  yielded  essen¬ 
tially  the  same  joint  first  order  analysis  and  model  parameters,  yet  the 
particular  time  structure  differences  of  the  two  sample  records  produced 
these  simulation  differences. 

We  note  also  in  these  three  simulations  that  the  discrepancies  be¬ 
tween  linearized  predictions  and  actual  test  results  are  reduced  for 
decreasing  estimator  bandwidths  (increasing  observation  intervals)  as 
we  would  expect.  Thus  the  quiet-night  conditions,  which  we  noted  in 
Chapte:  2  exhibited  a  long-time  correlation  in  the  A(t^)  estimates,  pro¬ 
duced  a  performance  difference  at  the  largest  estimator  bandwidth  which 
disappeared  for  the  other  two  bandwidths.  A  similar  comment  applies 
to  the  right-hand  (afternoon)  tropical  noise  record.  The  cross-over  in 
the  hard  limiting  versus  clipping  performance  for  frontal  noise,  pre¬ 
dicted  by  the  ZNL  I.  F.  curves  and  observed  in  the  scalar  simulations, 
does  not  occur  for  the  two -state  model  although  again  the  differences 
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decrease  with  decreasing  estimator  bandwidth. 

In  contrast  to  these  observations,  we  see  that  the  quiet  noise  condi¬ 
tions,  which  had  no  significant  multiple  discharge  activity,  produce 
results  agreeing  with  the  linearized  analysis  for  all  signal  models  and 
Q/SNE  conditions.  We  can  conclude  then  that  the  performance  prediction 
discrepancies  are  directly  attributable  to  this  noise  time  structure 
wherein  successive  samples  are  not  statistically  independent,  although 
uncorrelated. 

Two  representative  plots  of  the  time  evolution  of  the  estimate  error 
are  shown  in  Figure  4-9  where  the  relative  differences  in  performance 
given  by  the  error  statistics  can  be  seen. 

The  results  of  the  two -state  simulation  indicate  that  the  linearized 
ZNL  analysis  is  valid  for  most  noise/signal  conditions  of  interest,  although 
this  begins  to  break  down  for  short  estimator  bandwidth  relative  to  the 
noise  time  structure.  Using  (4. 43),  the  upper  left  simulation  has  approx¬ 
imately  50  samples  per  estimator  time  constant,  which  for  the  1  kHz 
sampling  rate  used  for  all  of  the  noise  records  corresponds  to  a  50-100 
msec  observation  interval.  We  recall  from  Chapters  2  and  3  that  mul¬ 
tiple  discharge  phenomena  have  typical  lengths  in  the  200-400  msec 
region.  Thus  significant  growth  of  signal  state  uncertainty  (refer  to 
section  4. 4.  2  discussion)  can  occur  during  such  a  noise  burst  which  is 
not  optimally  compensated  for  by  the  ZNL  estimators, or,  reflected  in  the 
linearized  analysis.  Corresponding  samples  per  time  constant  figures 
for  the  other  two  simulations  are  200  and  628,  indicating  that  the  effect 
of  multiple  discharges  becomes  progressively  less  significant  with  the 
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Figure  4-9.  Time  Evolution  of  Estimation  Error  for  Two-State  Signal 
Model 
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performance  results  noted.  The  general  improved  performance  of  the 
time-varying  estimator  again  validates  the  joint  channel  model  concept 
and  also  indicates  under  what  noise/bandwidth  conditions  the  use  of  this 
model  and  filter  results  in  significant  performance  improvement.  Gen¬ 
erally, this  improvement  is  significant  for  more  complex  signal  dynamics 
and/or  observation  intervals  that  are  comparable  to, or  less  than, the 
longest  interval  of  noise  statical  dependence. 

4.  7  Summary  and  Remarks 

4.  7. 1  Summary  of  Navigation  Receiver  Design  Results 

In  this  chapter  we  have  considered  the  problem  of  estimating  the 
phase  of  a  received  rf  carrier  when  the  observations  are  corrupted  by 
additive  atmospheric  noise.  The  phase  of  the  rf  signal  may  be  consid¬ 
ered  as  a  random  signal  process  representing  the  geographic  position 
of  the  receiver  relative  to  the  transmitting  station.  As  such,  this  prob¬ 
lem  is  the  most  fundamental  sub -problem  in  the  design  of  a  navigation 
receiver.  A  linearized  approximation  to  the  sampled  received  phase 
was  used,  and  three  types  of  estimator  structures  were  considered;  a 
linear  time -invariant  estimator  matched  to  the  received  noise  power 
and  signal  model,  a  time -invariant  estimator  using  a  zero  memory 
nonlinearity  (ZNL)  acting  on  the  loop  error  signal  to  attenuate  large 
noise  excursions,  and  a  near  optimum  time-varying  estimator  utilizing 
samples  of  the  noise  envelope  from  an  adjacent  frequency  channel  as  a 
measure  of  the  time -varying  power  level  of  the  received  noise. 

The  performance  of  the  ZNL  type  of  estimator  was  analyzed  by 
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replacing  the  ZNL  with  an  equivalent  gain  and  replacing  the  received 
noise  with  the  difference  between  the  ZNL  output  and  actual  error  signal. 
For  the  case  of  a  small  error  signal,  this  loop  could  be  optimized  through 
the  use  of  a  linear  time -invariant  estimator  following  the  ZNL,  which 
was  matched  to  the  effective  noise  power  at  the  ZNL  output  and  the  sig¬ 
nal  model.  This  approximate  analysis  allowed  complete  estimator  per¬ 
formance  to  be  specified  from  an  analysis  of  the  ZNL  based  on  the  fi  -st- 
order  pdf  of  the  noise  and  error  signal.  This  analysis  was  carried  out 
using  the  first-order  noise  model  from  Chapter  3.  Four  types  of  ZNL's 
were  analyzed  in  this  manner,  the  optimum  (MMSE)  ZNL,  hard  limiter, 
clipper  and  hole  puncher.  Based  on  this  analysis  two  types  of  ZNL's 
wei  3  considered  most  practical  for  actual  use,  a  hard  limiter  and  clipper. 

The  performance  of  the  third  type  of  estimator,  the  time-varying 
Kalman-Bucy  filter,  c^uld  not  be  quantitadvely  analyzed  due  to  the  non¬ 
linear  dynamic  covariance  equations  involved.  The  qualitative  operation 
of  the  loop  was  discussed  and  shown  to  provide  essentially  the  same 
behavior  as  the  ZNL  in  terms  of  attenuating  large  noise  excursions  and, 
in  addition,  providing  a  dynamic  compensation  for  the  loss  of  signal 
information  during  large  atmospheric  noise  bursts.  An  optimistic  lower 
bound  on  estimator  error  was  related  to  the  variance  of  the  background 
component  of  atmospheric  noise  which  is  of  constant  power  level  and 
uncorrelated  between  frequency  channels.  A  more  practical  lower  bound 
was  suggested  to  be  that  predicted  by  the  MMSE  ZNL  estimator,  although 
a  formal  proof  of  this  conjecture  was  not  obtained. 

Four  estimator  structures,  linear  time-invariant,  hard  limiting, 
clipping  and  me -varying  were  simulated  on  a  computer.  Actual 
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atmospheric  noise  sample  records  were  used  to  provide  the  noise  com¬ 
ponent  and  a  computer-generated  gaussian  random  process  was  used  to 
represent  the  signal.  The  results  of  these  simulations  showed  that  the 
linearized  ZNL  analysis,  based  on  the  atmospheric  noise  model,  pro¬ 
vided  quite  accurate  predictions  of  actual  performance.  Where  the  actual 
error  variance  was  found  to  be  greater  than  that  predicted,  the  optimum 
time-varying  estimator  generally  provided  improved  performance,  ap¬ 
proaching  the  optimum  predicted.  The  conditions  required  to  realize 
this  improvement  over  a  time -invariant  ZNL  estimator  were  complex 
signal  dynamics  (a  two-state  signal  model  as  compared  to  a  scalar)  and 
an  effective  estimator  time  constant,  or  observation  interval,  on  the 
same  order,  or  less  than.the  longest  time  interval  of  correlated  noise 
statistics,  approximately  0.  5-1  second.  These  application  results  verify 
our  hypothesis  from  Chapters  2  and  3  that  the  most  significant  charac¬ 
teristic  of  atmospheric  noise  is  the  first-order  pdf,  with  the  noise-burst 
time  structure  of  second-  or  third-order  significance.  Where  the  latter 
is,  or  thought  to  be,  important,  it  can  be  optimally  accounted  for  through 
the  use  of  the  joint  channel  noise  model  and  the  concept  that  the  noise 
waveform  can  be  modeled  as  a  time -varying  gaussian  process,  which 
leads  to  the  optimum  time -varying  estimator. 

4.  7.  2  Remarks  on  Re ceiver  System  Design 

Prior  to  offering  some  general  conclusions  on  the  suitability  and 
realization  of  the  three  basic  types  of  PLL's  we  have  studied,  we  shall 
consider  two  related  topics. 
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a)  Other  ZNL:s.  The  generally  excellent  results  we  have  obtained 
in  analyzing  three  specific  types  of  ZNL's,  operating  in  all  types 
of  atmospheric  noise  conditions,  suggests  that  other  ZNL  types, 
of  interest  in  a  particular  application,  can  be  similarly  analyzed 
with  confidence  in  the  resulting  performance  predictions.  The 
use  of  the  ZNL  Improvement  Factor  and  noise  model  parameter 
sets  provide  the  designer  the  tools  necessary  to  do  this. 

b)  RF  Bandwidth  Effects.  If  we  compare  the  effective  noise  reduc¬ 
tion  factor  from  the  plots  of  simulation  results  with  the  changes 
in  normalized  rms  noise  level  given  in  Table  3-1,  we  see  that 
the  best  ZNL  estimator  will  achieve  a  performance  nearly  equal 
to  that  caused  by  the  background  gaussian  noise  component  (except 
for  the  severe  and  atypical  frontal  noise  conditions).  In  other 
words,  the  estimator  will  be  nearly  insensitive  to  increases  in 
noise  level  caused  by  changing  weather/noise  conditions.  We 
recall  from  section  3.  9.  2  that  the  effect  of  reducing  the  bandwidth 
of  the  rf  bandlimiting  filter  is  to  increase  the  apparent  back¬ 
ground  gaussian  noise  component  for  increasingly  severe  noise 
conditions  ( px  increasing).  Thus,  the  ZNL  estimators  will  not 
maintain  constant  weather/noise  performance  with  rf  bandwidths 
significantly  less  than  1  kHz.  Rather  they  will  tend  towards  the 
degradation  in  performance  seen  for  an  ideal,  linear,  time- 
invariant  system.  While  our  experimental  investigations  have 
not  identified  the  rf  bandwidth  at  which  these  effects  become 
significant,  our  application  results,  based  on  the  noise  model, 
indicate  the  largest  effective  noise  reduction  one  can  obtain  with 
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an  optimal  bandwidth  selection,since, there  are  no  significant 
changes  in  either  model  or  application  results  for  bandwidths 
larger  than  1  kHz. 

4.  7.  2. 1  Linear  Receivers 

A  general  result  from  the  analysis  and  tests  of  navigation  PLL's  in 
this  chapter  is  that  one  should  not  employ  a  wide  dynamic  range,  linear 
PLL,  for  signals  corrupted  by  atmospheric  noise.  Where  the  other 
important  aspects  of  a  design  problem  (such  as  interference  suscepti- 
bility  *  )  dictate  the  use  of  a  linear  type  system,  some  sort  of  clipper 

nonlinearity  should  also  be  employed.  Our  results  show  that  this  pro¬ 
vides  near  optimum  operation  for  most  signal -to -noise  ratios,  signal 
models,  and  types  of  atmospheric  noise.  As  we  have  noted,  the  imple¬ 
mentation  of  the  adaptive  threshold  clipper  that  we  have  used  is  quite 
easy,  although  any  similar  type  of  threshold  setting  procedure  would 
probably  perform  nearly  as  well.  The  only  other  parameter  required 
in  such  a  system  to  achieve  optimum  performance  for  a  given  signal 
model  (or  alternatively,  in  a  classical  context  as  defined  in  section  4.2.3, 
to  maintain  a  fixed  loop  bandwidth)  is  the  effective  noise  power,  which 
can  be  estimated  from  the  variance  of  the  ZNL  output  sequence.  Thus 
an  intuitive  and  reasonable  modification  to  a  basic  linear  PLL  will  pro¬ 
vide  near  optimum  performance  for  atmospheric  noise. 

4.  7.  2.  2  Hard  Limiting  Receivers 

The  great  attraction  of  a  hard  limiter  ZNL,  as  we  have  employed 
it  here,  is  chat  the  PLL  can  be  almost  entirely  realized  with  digital 
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circuits.  A  stable  analog  oscillator,  low  level  preamplifier,  and  limiter 
circuits  are  the  only  analog  elements  required  in  the  general  case,  although 
other  variations  such  as  analog  feedback  to  a  voltage -controlled  oscillator 
can  be  employed.  The  hard  limiting  receiver  provides  near  optimum 
performance  in  most  cases,  especially  those  with  longer  observation 
intervals  which  are  the  cases  of  greatest  practical  interest.  Another 
advantage  that  is  often  associated  with  hard  limiting  PLL's  is  that  they 
do  not  appear  to  require  automatic  gain  control  loops  (AGC).  However, 
this  is  not  quite  the  case.  The  function  of  AGC  in  the  linear  receiver 
(with  or  without  clipper)  is  to  stabilize  the  effective  loop  error  gain  (K) 
at  a  prescribed  value  which  is  optimum  for  a  given  signal  model  and 
noise  power,  or  in  the  classical  design  sense,  results  in  a  prescribed 
loop  bandwidth.  The  hard  limit  PLL  requires  exactly  the  same  sort  of 
information  for  exactly  the  same  purpose,  namely,  the  effective  gain  c 
of  the  limiter  ZNL  and  the  effective  signal-to-noise  ratio  to  specify  the 
optimum  loop  gain,  K.  We  see  from  the  expressions  (4.  29)  and  (4.  30) 
defining  these  two  parameters  that  they  depend  only  upon  the  probability 
density  of  the  loop  error  process,  z(k+l),  evaluated  at  zero,  pz(0).  This 
can  be  estimated  in  a  practical  application  by  using  a  window  discrimin¬ 
ator  with  a  small  but  finite  width, centered  at  zero  volts,and  a  counter. 

This  would  determine  the  average  occupancy  of  this  increment,  by  the 
input  analog  waveform,  at  the  VCO  zero  crossing  time  (sample  time). 


$  45  48 

While  hard  limiters  tend  to  be  self-adaptive  over  a  limited  range  ’ 
this  variation  does  not  result  in  optimal  compensation  nor  does  it  con¬ 
sider  the  additional  problem  of  adaptation  to  the  changing  character  of 
of  atmospheric  noise. 
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Thus,  a  hard  limiting  system  can  be  designed  to  provide  near  optimum 
performance  for  most  atmospheric  noise  and  SNR  conditions  without 
sacrificing  realization  simplicity.  A  "choice"  between  a  hard  limit  or 
linear -with-clipping  type  PLL  for  the  navigation  application  would  seem, 
then,  to  depend  upon  factors  other  than  the  atmospheric  noise  perform¬ 
ance  that  we  have  investigated. 

4.  7.  2.  3  Optimum  Time -Varying  Receiver 

The  optimum  time -varying  PLL  that  we  have  considered  in  this  work 
is  definitely  more  complex  than  either  of  the  two  ZNL  types,  requiring 
both  the  co-channel  envelope  information  and  small  digital  computer  to 
compute  the  covariance  propagation  and  resulting  time-varying  gain, 
K(k+1).  The  latter  capability  is  often  found  in  modern  receiver  systems, 
especially  integrated  navigation  systems  that  combine  inertial  platforms 
or  Doppler  radar  systems  with  radio  aids.  The  co-channel  information 
can  be  provided  by  a  relatively  unsophisticated  algorithm  to  determine 
the  two  points  on  the  c r  i  curve  that  were  used  to  establish  E(k+1)  - 

yiv0 

R(k+1)  mapping.  Such  an  optimum  system  should  provide  excellent 
adaptibility  to  any  type  of  atmospheric  noise  condition  found  throughout 
the  world.  However,  there  is  very  little  performance  advantage  in  the 
system  for  estimator  time  constants  greater  than  several  seconds,  and 
these  are  the  ranges  of  greatest  practical  interest.  We  conclude,  there¬ 
fore,  that  the  principal  use  of  the  optimum  structure  based  on  the  joint 
channel  noise  model  has  been  to  demonstrate  the  degree  of  optimality  of 
other,  more  practical  systems,  over  the  range  of  signal  and  noise  con¬ 
ditions  that  we  have  investigated. 


Chapter  5 


SUMMARY  AND  SUGGESTIONS  FOR  FURTHER  WORK 

5. 1  Summary 

The  work  reported  here  has  involved  three  distinct  phases;  experi¬ 
mental  observations,  mathematical  modeling,  and  use  of  the  model  in 
a  practical  design  problem.  As  we  have  stressed  at  several  points, the 
end  use  of  the  model  is  an  engineering  tool,  and  this  has  provided  moti¬ 
vation  and  guidance  as  to  approximations  and  characteristics  that  were 
more,  or  less,  important  in  the  previous  efforts.  Three  significant 
results  of  the  experimental  and  modeling  efforts  can  be  summarized 
as  follows 

1)  The  model  provides  a  first  order  probability  density  for  a  band- 
limited  atmospheric  noise  waveform  or  noise  envelope,  and  a 
first  order  joint  probability  density  for  similar  waveforms 
observed  in  disjoint  frequency  channels. 

2)  The  time  structure  of  the  noise  model  approximates  at  least 
one  important  aspect  of  the  long  time  statistical  dependence  of 
atmospheric  noise  waveforms  (beyond  the  correlation  time  of  the 
bandlimiting  filter). 

3)  Comparison  of  the  noise  model  with  extensive  measurements 

of  bandlimited  atmospheric  noise  shows  that  the  model  matches 
important  noise  characteristics  over  a  broad  range  of  noise  con¬ 
ditions.  This  comparison  provides  a  table  of  model  parameters 
describing  this  range  of  noise  behavior. 
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The  use  of  the  atmospheric  noise  model  in  this  work  was  restricted 
to  a  specific  problem  and  considered  only  certain  sub-specifics  within 
that  problem,  such  as  the  type  of  zero  memory  non-linearity  to  use  in 
a  filter,  A  discussion  of  these  details  is  given  in  Chapter  4.  On  a 
broader  level  we  believe  there  are  two  significant  results  from  this 
application  of  the  noise  model: 

1)  A  useful  approximation  for  atmoshperic  noise  is  to  consider 
samples  spaced  beyond  the  bandlimiting  filter  correlation  time 
to  be  statistically  independent.  Under  this  assumption  the  first 
order  probability  density,  for  all  types  of  noise  conditions,  pro¬ 
vides  the  information  necessary  to  determine  the  performance 
of  various  types  of  receiver  structures.  This  independence 
assumption  and  use  of  the  first  order  model  yielded  the  ZNL 
filter  analysis  of  Chapter  4. 

2)  The  joint  channel  noise  model  and  the  concept  of  modeling  the 
noise  as  a  time -varying  Gaussian  process  allows  one  to  deter¬ 
mine  the  near  optimum  performance  of  a  receiver  operating  in 
all  types  of  atmospheric  noise.  Even  if  such  a  design  would  not 
be  considered  practical  for  implementation,  it  provides  a  funda¬ 
mental  limit  for  comparison  with  other  simplified  designs. 

5.  2  Suggestions  for  Additional  Work 
5.  2.  1  Noise  Model 

As  noted  at  the  end  of  Chapter  3  the  most  significant  questions  for 
further  experimental  research  on  this  type  of  noise  model  involve  fre¬ 
quency  dependent  effects.  The  most  important  would  be  a  determination 
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of  the  effect  of  rf  bandwidth  on  the  model  parameters.  While  our 
results  are  for  relaii  dly  broad  bandwidths  that  allow  one  to  deter¬ 
mine  the  maximum  performance  improvement  possible  with  atmo¬ 
spheric  noise,  a  designer  must  know  at  what  point  bandwidth  reduction 
begins  to  penalize  t.iese  improvements.  Other  significant  questions 
might  be  listed  as  follows  • 

1)  Explore  the  joint  channel  model  in  more  detail  at  VLF  where 
the  impulsive  return  stroke  is  more  important  than  the  leader 
noise  burst. 

2)  Determine  the  effectiveness  of  the  noise  model  at  frequencies 
above  LF.  Since  the  noise  at  these  frequencies  is  increasingly 
dominated  by  the  noise  burst  properties  of  the  discharge,  we 
would  expect  that  the  model  would  be  increasingly  effective. 

3)  Determine  the  effect  on  the  joint  channel  model  of  separation 
between  channels  at  VLF,  LF,  and  above. 

The  measurement  suggested  here  could  probably  be  most  conveniently 
made  by  an  instrumentation  system  using  frequency  shifting  techniques 
to  standardize  filter  shapes  at  the  intermediate  frequency.  We  avoided 
this  approach  in  our  work  because  of  the  difficulty  in  realizing  extremely 
wide  dynamic  range  multipliers  to  accomplish  the  shifting  operation.  Our 
results  show  that  the  application  of  the  noise  model  is  not  critically  sen¬ 
sitive  to  the  large  noise  amplitudes  and  hence  a  system  of  reduced  dynamic 
range  that  introduced  a  zero  memory  distortion  (no  memory  over  the 
sampling  interval)  at  these  amplitudes  would  not  seriously  compromise 
the  resulting  data  analysis  and  modeling. 


5.  2.  2  Application  to  Communication  Receivers 
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The  low  frequency  electromagnetic  spectrum  accommodates  a 
greater  number  of  digital  communication  systems  than  radio  naviga¬ 
tion  systems.  Most  of  these  digital  signalling  schemes  employ  bit 
lengths  of  the  order  of  20-50  msec.  This  time  period,  the  length  of 
the  observation  interval  in  which  the  receiver  must  decide  which  of 
m  signals  was  sent,  is  in  the  range  where  we  found  that  the  optimum 
waveform  estimation  receiver  offered  the  greatest  improvement  over 
simple  non-linear  processing  techniques.  We  suspect  that  similar  and 
perhaps  more  dramatic  improvements  are  possible  in  the  digital  re¬ 
ceiver  detection  problem. 

Just  as  in  the  estimation  application,  the  model  of  the  noise  pro¬ 
cess  as  a  time-varying  Gaussian  process  provides  immediate  application 
of  the  general  results  of  detection  theory  that  include  time -varying 
Gaussian  power  levels.  In  the  simplest  form  of  detection  problem, 
deciding  whether  a  known  signal  is  present  or  not,  we  can  draw  the 
form  of  the  optimum  detection  receiver  directly  from  the  insight  pro¬ 
vided  by  the  joint  channel  noise  model.  In  Figure  5-1  we  see  in  the 
upper  portion  a  basic  correlation  receiver  for  this  problem.  The 
received  waveform,  r(t),  is  correlated  against  a  replica  of  the  trans¬ 
mitted  signal  and  the  result  integrated  over  the  observation  interval  to 
form  the  test  statistic.  The  block  labeled  f(z)  is  a  linear  transfer 
function  in  the  linear  form  of  the  receiver,  which  implicitly  assumes 
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the  atmospheric  noise  is  white  and  Gaussian.  In  typical  atmospheric 

noise  applications  this  function  is  a  clipper  or  hole  puncher  or,  it  can 

be  specified  in  a  more  optimal  manner  using  the  results  of  various 
43  44  50 

studies  ’  ’  concerned  with  the  optimal  non-linearity  for  statis¬ 

tically  independent  samples  of  non-Gaussian  noise.  The  decision 

37 

threshold  is  then  determined  from  the  criterion  to  be  used  (Bayes, 
Neyman-Pearson,  Min-Max  etc)  and  the  effective  noise  and  signal 
power.  These  latter  parameters  must  be  calculated  in  a  manner  simi¬ 
lar  to  the  analysis  used  in  Chapter  4  for  the  ZNL  estimator.  One  would 
use  the  first  order  pdf  of  the  noise  (thus  using  an  ensemble  average 
of  noise  waveforms  for  that  noise/weather  condition)  and  compute  the 
signal  suppression  and  the  effective  noise  power  at  the  output  of  f(z). 

The  lower  portion  shows  the  nearly  optimal  use  of  the  pilot  chan¬ 
nel  noise  information.  The  pilot  channel  filter,  F2(f),  is  specified  to 
have  exactly  the  same  lowpass  envelope  response  as  the  signal  channel 
filter  Fj(f).  This  eliminates  the  problem  of  time  offset  found  in  our 
"worst -case"  noise  observation  of  Chapters  2  and  3  and  should  pro¬ 
vide  better  correlation  between  waveforms  than  our  experimental 


We  note  that  the  atmospheric  noise  component,  y(t),  has  a  correlation 
function  determined  by  the  bandlimiting  filter,  Fj(f).  However,  this  fil¬ 
ter  is  used  to  remove  adjacent  channel  interference  and  hence  r(t)  can¬ 
not  be  "whitened,"  the  y(t)  component  must  be  treated  as  white  noise  by 
the  receiver. 
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results.  The  non -linearitytg(v),  realizes  the  inverse  of  the  E(k+1)  — 

[R(k+l)j  '  mapping  of  Chapter  4  so  that  the  output  can  be  directly 

applied  to  the  indicated  multipliers.  This  results  in  an  effective 

received  signal,  r'(t),  that  has  a  nearly  constant  power  level  and  a 

nearly  Gaussian  noise  distribution  (see  section  3.4).  This  waveform 

is  then  correlated  against  a  new  s{t),  s'(t),  which  reflects  the  effect 

of  the  scaling  waveform  derived  from  the  pilot  channel.  The  result 

may  then  be  clipped  by  a  fixed  threshold  f(z)  to  include  the  fact  that 

the  joint  channel  noise  model  is  an  approximation  and  the  effective 

noise  power  component  of  r'(t)  is  not  exactly  Gaussian  (this  parallels 

the  time -varying  clipper  of  Chapter  4).  After  integration  across  the 

observation  interval  to  form  the  test  statistic,  the  decision  is  made 

2 

using  the  effective  signal  power,  s^ ,  received  during  that  interval  and 
the  nearly  constant  noise  power  of  r'(t). 

The  duality  of  the  estimation  problem  of  Chapter  4  and  the  detec¬ 
tion  scheme  proposed  here  is  evident.  The  time -invariant  ZNL  design 
and  analysis  uses  the  same  assumptions  and  type  of  first  order  pdf 
analysis.  The  optimal  time -varying  detection  receiver  parallels  the 
time -varying  estimation  receiver  in  the  following  ways;  the  inverse 
scaling  operation  attenuates  noise  excursions  in  the  same  way  as  the 
estimator's  use  of  R_1(k+1)  information  (section  4.4.  2)  and  the  compu¬ 
tation  of  effective  signal  power  for  each  observation  interval  optimizes 
the  decision  threshold  for  that  interval  just  as  the  time -varying  estimator 
optimally  weighted  each  sample  via  the  S(kH  |k+l)  term  of  the  gain 
matrix  K(k+1). 
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We  suspect  that  the  optimum  detection  receiver  will  provide  greater 
improvements  than  the  optimum  estimator  tests  indicated  because  of 
the  short  observations  intervals,  between  which  the  atmospheric  noise 
waveform  is  very  non-stationary,  and  the  fact  that  the  estimator  test 
records  were  based  on  samples  spaced  1  msec  apart,  whereas  the 
detection  receiver  can  make  effective  use  of  much  higher  information 
rates.  We  also  note  that  the  optimum  detection  receiver  employs  principally 
analog  multipliers  which  can  be  easily  realized  to  operate  at  these  VLF 
and  LF  frequencies.  Extension  of  the  suggested  detector  to  actual 
FSK  signalling  schemes  should  be  straightforward.  We  can  also  en¬ 
vision  further  uses  of  the  joint  channel  information,  for  example,  to 
M*Iag"  bit  decisions  which  have  a  high  probability  of  error  based  on  the 
actual  signal -to -noise  ratio  for  that  observation  interval.  These  flags 
could  be  employed  in  following  decoders  for  error  correction  or  retrans¬ 
mission  requests. 
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Appendix  A 

INSTRUMENTATION  SYSTEM 


A..  I  Introduction 

This  appendix  provides  a  general  functional  description  of  the  instru¬ 
mentation  system  used  to  collect  information  on  atmospheric  noise  for 
this  research.  Detailed  circuit  schematics  of  the  equipment  are  not 
provided,  although  measured  performance  characteristics  of  certain 
critical  circuits  are  given. 

The  general  functions  of  the  system  can  be  categorized  as  follows: 

i)  Receive  and  process  low-frequency  atmospheric  noise  wave¬ 
forms  in  the  frequency  region  14  kHz  to  100  kHz. 

ii)  Sample  and  digitize  the  processed  waveforms  according  to  a 
prescribed  timing  sequence. 

iii)  Record  the  resulting  digital  data  in  a  manner  consistent  with 
computer  input  facilities  available. 

In  meeting  the  above  functional  requirements,  several  performance 
categories  were  important.  The  first  of  these  was  to  provide  a  dynamic 
range  in  the  analog  circuits  consistent  with  that  of  the  anticipated  atmo¬ 
spheric  noise  waveforms.  This  range  has  often  been  observed  in  excess 
of  100  dB;  however,  this  equipment,  due  to  restrictions  in  the  practically 
achievable  dynamic  range  of  high  speed  analog-to-digital  converters  (A/D) 
is  designed  for  80  dB  in  nonlinear  processing  circuits  and  100  dB  in  pre¬ 
ceding  linear  circuits.  A  second  major  consideration  in  the  analog  cir¬ 
cuits  was  to  insure  that  the  atmospheric  noise-waveforms  were  not 
obscured  by  man-made  interference,  either  communication  signals  trans¬ 
mitted  throughout  the  frequency  region  of  interest  or  noise  generated 
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within  the  equipment.  The  latter  may  be  either  random  noise  in  the  low- 
level  analog  circuits  or  digital  switching  noise  generated  in  the  timing 
circuits,  and  hence  synchronous  with  the  sampling  patterns. 

Within  the  digital  portions  of  the  equipment  the  principal  performance 
requirement  was  the  generation  of  a  wide  variety  of  sampling/data  hand¬ 
ling  patterns,  varying  from  a  maximum  rate  of  100  kHz  to  a  minimum 
rate  of  0. 1  Hz.  A  second  major  requirement  was  the  capability  to  inter¬ 
face  this  wide  variation  in  data  acquisition  rate  with  the  data  recording 
device. 

A  final  major  performance  specification  was  the  inclusion  of  an 
operator  interface  system  that  provided  rapid  and  flexible  adjustment 
of  system  parameters  and  the  accurate  monitoring  and  display  of  vari¬ 
ous  signal  levels  so  that  the  available  system  capabilities  could  be 
individually  optimized  for  each  data  sequence.  This  overall  performance 
is  especially  important  when  we  recognize  that  the  atmospheric  noise  is 
stationary  over  a  period  often  measured  in  hours  and  that  a  wide  variety 
of  sample  records  had  to  be  collected  within  that  time  frame  to  charac¬ 
terize  the  noise  under  those  particular  conditions. 

An  overall  functional  diagram  of  the  equipment  is  shown  in  Fig¬ 
ure  A-l.  A  sub -functional  description  with  associated  performance 
measurements  is  given  in  sections  A.  2  and  A.  3  for  the  analog  and  digital 
portions,  respectively. 

A.  2  Analog  Signal  Processing 

Low-frequency  atmospheric  noise  fields  were  received  on  a  one- 
meter  loop  antenna  manufactured  by  Aerospace  Research  Inc.  A  single- 
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A-l.  Block  Diagram  of  Instrumentation  System 
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turn  primary,  coupled  to  a  seven-turn,  secondary  transformer,  drove 
a  low- input  impedance  preamplifier  via  a  series  resonating  capacitor. 
This  circuit  provided  a  single  pole  response  of  approximately  60  kHz 
bandwidth  that  was  tuned  to  either  100  kHz  or  36  kHz.  The  antenna  pre¬ 
amplifier  had  an  equivalent  noise  voltage,  referred  to  the  input,  of 
8.  5  nV/sl Hz.  This  level  was  approximately  6-10  dB  below  the  back¬ 
ground  atmospheric  noise  component,  at  100  kHz,  as  determined  by 
visual  observation  of  the  respective  waveforms  on  an  oscilloscope.  The 
preamplifier  output  was  matched  to  a  100-ohm,  balanced,  transmission 
line  which  also  supplied  amplifier  power. 

The  signal  was  then  coupled  to  an  array  of  L  -C  notch  filters  to  re¬ 
move  large  amplitude  communication  signals  present  throughout  the  low- 
frequency  bands.  These  filters  were  provided  by  the  U.S.  Coast  Guard 
Electronics  Engineering  Center,  Wildwood,  N.  J. ,  and  exhibited  maxi¬ 
mum  notch  depths  of  60  dB  with  notch  widths  at  the  -3  dB  points  of  2  kHz. 
The  notch  filter  output  was  coupled,  via  a  balanced  line,  to  a  differential 
input  FET  amplifier  followed  by  a  broadband  operational  amplifier  and 
bootstrapped  twin-tee  notch  filter  to  remove  a  strong  interfering  signal 
at  18  kHz,  which  was  outside  the  tuning  range  of  the  L-C  notch  filters. 
The  composite  performance  of  the  analog  system  to  this  point  is  shown 
in  Figure  A-2. 

The  received  noise  signal,  with  man-made  interference  reduced  to 
manageable  levels, was  then  coupled  to  several  types  of  bandpass  filters 
to  establish  the  various  noise  frequency  channels  of  interest.  At  the 
outset  of  the  research,  two  single-stage  bandpass  filters  of  a  balanced 
Butterworth  type  were  used  for  this  function.  The  center  frequency  and 
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bandwidth  of  these  filters  could  be  conveniently  changed  with  a  plug-in 
resistor -capacitor  network.  The  inductors  used  were  the  same  high-Q 
pot  core  types  used  in  the  notch  filters,  and  they  permitted  a  maximum 
filter  Q  of  approximately  500.  As  the  research  progressed  and  certain 
frequency  channels  were  selected  as  providing  representative  data,  these 
filters  were  specialized  to  fix-tuned,  multi-stage  types  of  the  Tchebycheff 
or  Lerner  characteristic.  Figure  A-3  shows  the  frequency  response  of 
the  various  filters. 

The  bandlimited  noise  waveforms  could  be  connected  to  any  combi¬ 
nation  of  four  processing  channels;  two  wideband  amplifier/ samplers, 
or  two  linear  rectifiers  followed  by  lowpass  filters  or  integrate/ hold/ 
clear  circuits.  The  latter  provided  wideband  noise  envelopes  or  the 
A(t)  sample  series  discussed  in  Chapter  2.  The  wideband  samplers 
were  conventional  diode  ring  bridges  constructed  from  matched  quads 
of  Schottky  barrier  diodes.  The  ring  was  driven  by  two  dc  coupled, 
pulsed  current  sources  of  600nsec  nominal  width,  included  a  current 
balance  adjustment  and  a  variable  delay  adjustment.  The  sampling  gate 
analog  input  was  driven  by  a  low-impedance,  wideband,  operational 
amplifier  and  had  a  net  slew  rate  at  the  holding  capacitor  of  24  volts/ 
nsec.  This  rate,  coupled  with  the  fact  that  the  samples  were  taken  every 
10  nsec  and  the  noise  waveform  was  correlated  over  at  least  50  nsec, 
guaranteed  tracking  of  the  samples  over  the  full  ±10  volt  dynamic  range 
of  the  analog  system.  The  two  diode  bridge  current  adjustments  allowed 
independent  balancing  of  the  "on"  currents  and  the  turn-off  transient 
caused  by  both  sources  not  turning  off  at  exactly  the  sa  - e  time.  This 
balancing  reduced  absolute  offset  errors  in  the  resulting  sample  to  1  mV. 


Figure  A-3. 
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The  linear  rectifiers  utilized  germanium  diodes  in  the  feedback  loop 
of  an  operational  amplifier  to  provide  a  linear  range  of  7b  dB  when  prop¬ 
erly  trimmed.  The  rectifier  output  was  connected  to  a  two- stage,  3  dB 
ripple,  Tchebycheff  lowpass  filter  with  20  kHz  corner  frequency  to  recover 
the  envelope  from  the  rf  noise  waveform.  This  was  fed  to  the  wideband 
sampler  when  envelope  sample  records  were  required.  The  rectifier 
output  was  also  connected  to  a  gated  integrator  to  generate  the  A(t.) 
sample.  The  input  to  the  integrator  was  disconnected  by  an  FET  switch 
at  the  end  of  the  integration  period  and  the  output  held  for  25  p.sec  for 
A/D  conversion,  followed  by  a  75-jisec  clear  period  when  the  feedback 
capacitor  was  discharged  by  a  second  FET  switch.  This  "service  period" 
was  a  fixed  100  psec  for  any  specified  integration  period. 

The  four  possible  sampled  signal  sources  could  then  be  connected 
to  the  A/D  converter  via  four  FET  switches  and  a  wideband  buffer  ampli¬ 
fier  with  300 -nsec  settle  time.  These  switches  could  be  programmed 
to  scan  the  signal  sources  in  any  periodic  set  of  1,  2,  or  4  sources.  The 
switch  states  were  decoded  from  digital  control  sequences  and  the 
switches  were  activated  at  the  beginning  of  a  10-psec  sample/digitize/ 
store  time  block  so  that  transients  had  decayed  prior  to  initiation  of 
A/D  conversion.  The  A/D  converter  was  an  Analogic  AN-271.5M,  15-bit, 
successive  approximation  type  with  total  conversion  time  of  8.  5  psec  in 
our  application.  The  digital  output  was  in  parallel  l's  complement  and 
was  converted  to  sign  and  magnitude  prior  to  storage. 

The  entire  analog  system  was  connected  via  an  FET  signal  multiplex 
array  to  a  monitor  system.  Any  signal  point  (e.  g. ,  sampler  input, 
sampler  output,  rectifier  output,  etc.)  could  be  selected  for  monitoring. 
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resulting  in  signal  availability  at  a  front  panel  jack  for  oscilloscope  dis¬ 
play  and  processing  by  a  metering  circuit  for  average  value  of  the  signal 
magnitude  and  peak  signal  value.  These  parameters  were  displayed 
via  a  front  panel  digital  voltmeter.  The  peak  value  circuit  used  an  FET 
input  threshold  detector  that  provided  acquisition  of  a  lO-jiusec,  12-volt 
pulse  to  5%  accuracy,  with  a  hold  time -constant  of  10  minutes.  The 
combination  of  the  average  and  peak  value  circuit  allowed  the  operator 
to  monitor  these  two  parameters  of  the  wide  dynamic  range  atmospheric 
noise  waveforms  to  optimize  the  use  of  the  A/D  converter  dynamic  range. 
The  meter  system  was  also  used  for  calibration  of  dc  offsets  through¬ 
out  the  analog  processing  system.  In  addition  to  this  analog  monitoring, 
front  panel  lights  displayed  a  2-second  pulse  whenever  any  of  the  analog 
circuits  following  the  L-C  notch  filters  overloaded.  Finally,  the  A/D 
digital  output  was  displayed  in  sign  and  magnitude  format  for  real  time 
feedback  to  the  operator  of  the  effective  signal  range  and  for  use  in  cali¬ 
bration  of  converter  offset  and  gain. 

Figure  A-l  shows  the  various  timing  signals  provided  the  analog 
system  from  the  digital  system  and  also  an  analog -digital  path  labeled 
"level  discriminate."  This  was  used  for  conditional  sampling  patterns 
as  described  in  Chapter  2  and  consisted  of  two  discriminators.  The 
first  was  a  threshold  detector  connected  to  the  integrator  output.  This 
supplied  a  state  change  whenever  the  integrator  output  exceeded  a  pro¬ 
grammable  analog  threshold.  The  second  line  provided  a  digital  output 
from  a  window  discriminator  which  changed  state  only  when  the  sampled 
value  from  one  of  the  wideband  samplers  fell  within  a  programmable 
analog  voltage  window.  This  latter  circuit  had  a  minimum  window 


resolution  of  10  mV. 
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The  analog  signal  processing  system  utilized  state-of-the-art  de¬ 
vices  in  conventional  circuitry  that  was  carefully  optimized  to  provide 
the  required  dynamic  range,  operator  feedback  and  flexibility  to  con¬ 
duct  an  evolving  experimental  program  in  which  goals  and  techniques 
were  expected  to  change  as  information  and  insight  developed. 

A.  3  Digital  Timing  and  Control 

The  digital  timing  and  control  system  was  organized  as  a  synchro¬ 
nous  machine  on  a  basic  period  of  1C  nsec,  counted  down  from  the  1-MHz 
crystal  oscillator.  The  first  3  nsec  of  this  period  were  devoted  to  sam¬ 
pling,  signal  multiplexing  and  the  start  of  A/D  conversion.  The  4-5  nsec 
time  period  was  reserved  for  reading  data,  stored  in  core,  for  output  to 
the  paper  tape  punch.  The  last  2  nsec  were  reserved  for  formatting 
data  from  the  converter,  temporary  storage  and  writing  into  core  mem¬ 
ory.  Sampling  patterns  were  generated  as  multiples  of  this  10 -nsec 
period.  Sampling  could  be  specified  as  continuous  with  spacing,  T  ,  of 
10,20,50,  100,200  nsec,  etc.,  to  10  seconds,  or,  in  groups  of  size  2n. 
Group  spacing  was  independently  specified,  as  was  interior  sample 
spacing.  The  total  sample  record  could  be  specified  as  256  •  2n,  0  ^ 
n  9,  or  indefinite.  All  timing  functions  could  be  linearly  scaled  by  a 
factor  0.1  ^  a  «  1  for  synchronous  operation  at  other  band  center  fre¬ 
quencies.  This  was  accomplished  by  varying  the  input  oscillator  fre¬ 
quency  from  1  MHz  to  100  kHz. 

The  digital  process  control  was  a  dc  coupled,  finite  state,  machine 
which  could  be  placed  in  various  test/calibrate  modes  or  initialize/run 


-1*2  - 

mode.  In  the  latter  all  counters  and  check  circuits  were  cleared  and 
armed  by  front  panel  control  and  a  data  run  begun  by  pressing  the  "GO" 
button.  During  the  data  run,  the  state  of  various  digital  functions  was 
monitored  and  sampling  inhibited  in  the  event  of  the  machine  trying  to 
enter  a  disallowed  state  such  as  writing  new  data  over  data  that  had  not 
been  read  out  from  memory,  or  exceeding  the  programmed  number  of 
total  samples.  The  state  of  the  process  control  was  displayed  via  front 
panel  lights,  as  was  core  storage  remaining.  Synchronous  interrupt 
circuits  were  also  provided  to  allow  the  operator  to  stop  sampling  or  data 
outputting,  during  the  run,  without  requiring  re-initialization. 

A.  4  Data  Handling 

The  A/D  digital  output,  in  l's  complement,  was  converted  to  sign 
and  magnitude  and  stored  in  core  memory  in  one  of  three  forms.  The 
first  was  a  direct  transfer  of  the  entire  15 -bit  word  into  one  16 -bit 
core  memory  word.  This  mode  generally  provided  greater  sample 
resolution  than  necessary,  at  the  expense  of  data  throughput  to  the  re¬ 
corder,  and  was  only  used  for  applications  where  such  resolution  or  the 
sample  sign  was  required,  such  as  measuring  the  autocorrelation  of  the 
bandlimited  rf  noise  waveform.  The  second  mode  retained  only  the 
most  significant  8  bits  of  the  sample  magnitude,  packed  two  samples  per 
core  word,  and  was  used  principally  for  the  A(L)  sample  record  which 
was  unipolar  and  of  reduced  dynamic  raoge.  The  third  and  most  often 
used  mode  was  termed  "compressed"  and  retained  either  the  most  sig¬ 
nificant  7  bits  of  the  sample  magnitude  or  bits  6  through  12  if  bits  1 
through  5  were  all  zero  (where  bit  1,  Bj  -  MSB).  Bit^  of  the  stored 
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sample  indicated  whether  Bj-B^  or  Bg-Bjg  was  present.  This  mode 
of  operation  provided  varying  resolution  but  increased  dynamic  range 
for  maximum  data  throughput  and,  as  will  be  seen  in  Appendix  B,  was 
compatible  with  the  probability  density  estimation  algorithm  that  was 
used.  We  emphasize  that  while  the  second  two  modes  of  operation  sacri¬ 
ficed  digital  resolution  in  storing  and  transmitting  the  sample,  the  con¬ 
verter  always  operated  at  15  bit  resolution,  and  hence  the  precision  of 

each  7  or  8  bit  sample  was  that  of  the  original  15  bit  conversion.  Thus, 
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for  example,  the  resolution  of  mode  two  was  a  relative  2  but  this  bit 
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was  set  to  a  one  or  zero  with  a  precision  of  2  in  the  conversion  pro¬ 
cess. 

An  automatic  data  handling  test  system  was  built  in  and  exercised 
prior  to  each  major  data  recording  effort.  This  system  sent  a  staircase 
waveform,  of  4  steps,  to  the  A/D  converter,  with  binary  weighting  between 
steps.  This  was  converted,  stored  and  outputted  to  the  punch  in  all  data 
handling  modes  and  provided  easily  verifiable  bit  patterns  on  the  paper 
tape.  Since  the  waveform  was  derived  from  the  precision  reference 
source  of  the  A/D  converter,  it  tested  relative  converter  performance 
and  the  entire  data  format  and  processing  operation  through  to  the  punch. 

A. 5  Data  Storage  and  Output 

The  core  memory,  used  as  a  buffer  between  the  variable  sampling 
rates  and  paper  punch,  was  a  commercial  grade,  1. 5-^sec  access  time, 
memory  manufactured  by  Cambridge  Memories.  The  punch  was  a  Tally 
Corporation  Model  P-120  with  120,  8  bit,  character/second  capacity. 

The  punch  was  operated  on  a  "handshake"  basis  with  the  digital  process 
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control.  Punch  line  power  was  supplied  via  an  isolation  transformer 
and  the  data  and  logic  control  signals  were  transferred  via  an  optical 
isolator.  Thus,  there  was  no  common  ground  or  conducting  connection 
between  the  analog-digital  instrumentation  system  and  punch,  preventing 
ground  loop  conduction  ol  the  large  current  transients  generated  by  the 
punch  solenoid  drivers. 

24 

A  detailed  description  of  the  digital  system  is  given  by  Lee.  As 
we  noted  before  the  entire  instrumentation  system  was  characterized 
by  conservative  design  and  adaptability,  both  features  required  for 
study  of  noise  waveforms  whose  characteristics  were  not  known  in 
advance.  Each  subsystem  was  tested  during  assembly,  as  were  the 
analysis  programs  given  in  Appendix  B.  Finally,  a  one-month  period 
was  devoted  to  testing  the  entire  hardware/software  system  as  out¬ 
lined  in  Appendix  C. 


Appendix  13 
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DATA  ANALYSIS  ALGORITHMS 

Appendix  A  described  three  modes  in  which  the  binary  word,  repre¬ 
senting  the  data  sample,  were  recorded  on  paper  tape.  This  tape 
record  was  read  into  the  PDF- IX  computer  and  stored  on  magnetic 
tape  for  all  further  analysis  and  use.  The  tape  reading  subroutine 
unpacked  the  16  bit  word  and  rescaled  the  "compressed"  format  sample 
so  that  the  analysis  programs  treated  each  sample  as  a  14  bit  integer 
plus  sign,  with  only  the  value  of  the  least  significant  bit,  entered  at  the 
time  of  analysis,  indicating  the  actual  format  used  in  recording. 

B-l  Probability  Density  Analysis 

The  basic  frequency-of-occurrence  analysis  one  can  perform  on  a 

14  bit  integer  sample  record  (neglecting  sign)  is  to  count  the  number 

14 

of  occurrences  of  each  discrete  value,  or  which  there  are  2  -1  pos¬ 

sible  values.  In  terms  of  the  continuous  probability  density  of  the 

noise  waveform,  this  can  be  thought  of  as  sorting  all  samples  into 

- 1 4 

"bins"  of  normalized  width  2  .  The  probability  density  estimate 

based  on  each  bin  population  is  then 

Number  of  Samples  in  Bin  x. 

pd(x  )  =  — — 

1  (Total  Samples)  •  (LSB  Value) 

For  our  purposes  in  analyzing  atmospheric  noise  sample  records, 
which  have  a  large  dynamic  range  and  low  probability  of  occurrence 
in  the  largo  amplitude  region,  such  a  technique  leads  to  a  large 
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variance  in  the  individual  bin  estimates  and  is  prohibitive  in  terms  of 
memory  requirements  for  a  14  bit  maximum  word  length.  An  alter¬ 
nate  technique  that  smooths  the  probability  density  estimate  and  re¬ 
sults  in  a  small  memory  requirement  is  to  represent  each  sample  word 
as  a  binary  floating  point  number,N  =  a  •  2°,  where  'a"  is  a  4  bit 
number  and  0  <  b  <  10.  The  occurrence  of  each  sample  value  will  be 
represented  by  a  point  in  an  8  X  12  matrix  where  the  first  index  specifies 
the  8  possible  values  of  ''a"  (there  are  not  16,  except  at  b  =  0,  since 
a  given  value  o*  b  specifies  that  Bjj  ^  =  a^  =  1)  and  the  second,  the 
value  of  b.  Sorting  samples  in  this  manner  cr.n  be  easily  programmed 
with  a  rotate  instruction  to  fine*  the  leading  1  bit  (the  b  index)  and  then 
strobing  out  the  next  three  bits  for  the  a  index. 

The  principal  advantage  of  this  floating  point  analysis  algorithm 
is  that  it  provides  a  bin  width  which  expands  cxponentially.in  discrete 
steps, for  increasing  signal  amplitude.  viewed  another  way,  we  can 
say  that  it  provides  a  non -station ary  smoothing  function  for  the  proba  ¬ 
bility  estimate  since  it  averages  across  a  larger  number  of  the  basic 
bins  for  increasing  signal  amplitude.  The  probability  density  estimate 
for  this  procedure  is  given  by 

Nr  Samples  in  Bin  x. 

pd(x.)/volt  = - - - T". 

(Total  Samples)  •  (LSB  Value)  •  (2  ) 

where  LSB  value  is  the  value,  in  volts,  of  the  least  significant  bit  of 
the  data  format  used  on  the  record  being  analyzed.  For  use  in  cal¬ 
culating  sample  moments  the  signal  value  at  bin  center  is  used  (the 
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average  value  for  sign  and  magnitude  representation) 

x.  =  (a+l/2)  •  LSB  Value  •  2b. 

B-2  Probability  Density  Estimation  Error 
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It  is  well  known  that  the  probability  of  a  sample  occurring  in  a 
given  "bin"  is  a  binomial  trial  with  probability  of  occurrence 

r*x  +B 

p  =  \  °  W  p  (r\)  dtp  B  =  Bin  Width, 

*x 

o 


where  px(  • )  is  the  true,  continuous  probability  density  of  x„  With  the 
estimate  of  p  given  by  (Number  in  bin)/(Total  Samples  =  N)  the  vari  ¬ 
ance  of  the  estimate  is  (assuming  independent  samples) 

2 

2  PQ  _P  2  Pi 

“  n  ~  N  '  °norm  -  2  ”  pN  ’ 
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If  we  describe  px(  • )  as 
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X  <  X  , 
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and  consider  the  two  cases  of  constant  bin  width  and  the  exponentially 
increasing  bin  width  (where  we  approximate  the  discrete  increases  as 
a  continuous  relation) 
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This  qualitative  behavior  is  summarized  in  P'igure  B-l  where  we  see 
that  the  exponential  bin  width  associated  with  the  floating  point  number 
representation  provides  a  minimum  estimation  error  in  the  most  critical 
region  of  the  probability  density.  Greatest  weight  was  given  to  fitting 
the  model  to  the  estimated  density  in  this  region  (Chapter  3). 


Figure  B-l.  Probability  Density  Estimate  Error. 


B-3  Probability  Density  Surface  Analysis 


-199- 


The  probability  density  surface  analysis,  used  for  the  joint  chan¬ 
nel  data  records,  employed  a  floating  point  representation  for  the 
envelope  channel  sample  with  a  =  1  and  for  the  rf  channel  sample  with 
a  =  3.  This  resulted  in  a  4  X  13  X  13  matrix  to  represent  all  possible 
occurrences.  Density  estimation  and  moments  were  computed  as  in 
B-l  above. 

B-4  Correlation  Analysis 

All  sample  records  that  were  intended  for  correlation  analysis 
were  taken  in  continuous  bursts  where  the  burst  length  was  much 
longer  than  the  correlation  interval  of  interest.  The  basic  correlation 
technique  was  then  to  use  a  conventional  lagged  products  algorithm  on 
each  burst  and  then  average  each  correlation  point  across  all  bursts. 
Samples  from  different  signal  sources  were  interleaved  on  the  record 
without  synchronization  codes,  a  check  of  correct  data  transfer  was 
provided  by  total  sample  count  which  was  accurately  controlled  in  the 
instrumentation  system. 

If  we  define  the  following  quantities; 

A(i),  B(i)  -  Sample  Records, 

n  -  correlation  point  of  p^g(n),  N  -  Burst  Record  Length, 

K  -  maximum  value  of  n,  M  -  number  of  bursts, 


then  the  correlation  coefficient  was  estimated  as 
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pAB(n)  = 


A(i)  B{i+n) 

VA?'(i)  '  B2(i) 
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This  provides  an  unbiased  estimate  of  PAB(n),  however,  little  can  be 
said  about  estimate  variance  since  the  statistics  are  all  non-Gaussian. 
In  practice  the  entire  80000  point  sample  record  was  used  to  obtain 
the  best  stability  possible. 


Appendix  C 
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INSTRU MENTATION  TESTS 
C-l  Probability  Density  Tests 

The  first  test  that  was  conducted  of  the  instrumentation  hard¬ 
ware  and  analysis  software  was  a  uniform  distribution  test.  An 
asynchronous  square  wave  was  integrated  to  form  a  triangular  wave, 
which  when  sampled,  gave  rise  to  a  uniform  probability  distribution. 

The  principal  application  of  this  was  to  calibrate  slight  inequalities 
in  the  widths  of  the  A/D  converter  quantization.  These  were  apparent 
at  the  low  amplitude  estimates  of  the  probability  density  and  appeared 
as  a  noise  of  much  greater  variance  than  predicted  with  a  binomial 
assumption  (see  Appendix  B).  A  corrective  scale  factor  that  reduced 
the  uniform  probability  density  estimate  to  the  correct  value  was  then 
inserted  in  the  analysis  program  and  applied  to  all  data  analysis. 

The  second  test  involved  use  of  a  General  Radio  Random  Noise 
Generator  whose  output,  with  flat  power  spectral  density  across  the 
20  kHz  to  500  kHz  band,  was  applied  to  the  instrumentation  system  by 
transformer  coupling  to  the  loop  antenna.  The  resulting  bandlimited  noise 
waveform  was  sampled  at  10  times  the  filter  correlation  time  and 
the  resulting  probability  density  analysis  converted  to  a  cumulative 
distribution.  This  is  shown  below,  plotted  on  Gaussian  probability  paper 
where  a  Gaussian  cumulative  distribution  plots  as  a  straight  line. 
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The  last  probability  density  test  to  be  applied  consisted  of  an 
asynchronous  sine  wave  which  was  modulated  between  two  levels  of 

16  V  Peak-Peak  and 
200  mV  Peak-Peak.  This 
modulation  was  at  a 
100  Hz  rate  which  was 
asynchronous  with  both 
the  rf  sine  wave  and  the 
sampling  time-base.  This 
type  of  waveform  approxi¬ 
mated,  in  dynamic  range, 
the  type  of  behavior  anti¬ 
cipated  in  atmospheric 
noise  waveforms.  The  modulated  sine  wave  was  sampled  using  the 
"compressed"  data  mode  and  the  resulting  probability  density  estimate 

matched  quite  well  the  sum  of  two  l/Vl  -  xz/x2  functions  which  de- 

P 

scribe  such  a  switched , randomly  sampled, sine  wave. 

C-2  Correlation  Tests 

A  single  effective  test  of  the  A(t.)  estimate  generators  and  correla¬ 
tion  software  was  made  by  using  a  switched  attenuator,  as  described 
above,  to  modulate  the  power  level  of  a  white  Gaussian  noise  process. 
The  white  noise  was  then  bandlimited  at  two  different  center  frequencies 
and  the  A(t^)  estimates  of  the  short-time  power  level  were  formed  using 
1  msec  integration  times.  The  two  resulting  integrator  waveforms  are 


Figure  C-l.  Estimate  of  Gaussian  Cumula 
tive  Distribution 
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shown  in  Figure  C-2.  We  note  the  obviously  greater  variation  in  the 
estimates,  represented  by  the  peak  of  the  ramp,  in  the  narrowband 
waveform  compared  to  the  broader  bandwidth.  The  auto  and  cross¬ 
correlation  analysis  for  a  10,000  point  sample  record  is  shown  in  the 
same  figure.  The  waveform  variance  is  indicated  by  the  shaded  areas 
at  the  origin,  where  we  see  that  the  .75  kHz  bandwidth,  slightly  less  than 
used  in  the  final  data  recordings,  yields  a  normalized  variance  of  the 
power  level  estimate  of  0.  1.  In  comparison,  the  analysis  of  the  1  msec 
integration  time  data  indicated  variances  of  0.7-0.  85  and  we  would  con¬ 
clude  that  this  is  a  reasonably  accurate  measure  of  atmospheric  noise 
characteristics  rather  than  an  effect  caused  by  the  background  Gaussian 
noise.  Other  tests  verified  the  correlation  software  for  bipolar  data 
and  the  various  data  recording  modes. 


Figure  C-2.  Correlation  Test  Result 
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